REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  includirtg  the  time  for  reviewing  instructions,  searching  ex^tirig  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  Information.  Send  comments  reprding  this  burden  ^ 

rallection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate 

Davis  Highway,  Suite  1204,  Arlington^A  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Protect  10704-0188),  Washington,  DC  20503. 

1.  AGEWCY  USE  ONLY  (Leave  blank)  I  2.  REPORT  DATE  ^  |  3.  REPORT  TYPE  AND  DATES  COVERED 


28. Sep. 00 


TITLE  AND  SUBTITLE 


DISSERTATION 

I  5.  FUNDING  NUMBERS 


NONLINEAR  ULTRASOUND  IMAGE  MODELING:  DEVELOPMENT  OF  A 
COMPLETE  END-TO-END  MODEL  FOR  QUALITATIVE  AND  QUANTITATIVE 

ANALYSIS  OF  MEDICAL  ULTRASOUND  _ 

6.  AUTHOR(S)  ^  ~ 

MAJ  AYER  KEVIN  W 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

ROCHESTER  INSTITUTE  OF  TECHNOLOGY 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

CY00373 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
THE  DEPARTMENT  OF  THE  AIR  FORCE 
AFIT/CIA,  BLDG  125 
2950  P  STREET 
WPAFB  OH  45433 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Unlimited  distribution 

In  Accordance  With  AFI  35-205/AFIT  Sup  1 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 


14.  SUBJECT  TERMS 


15.  NUMBER  OF  PAGES 

123 

16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 


19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF 
OF  ABSTRACT  ABSTRACT 


ixnc  qualut  imemmsi  4 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


Nonlinear  Ultrasound  Inaage  Modeling:  Development  of  a  Complete  End-to-End  Model  for 
Qualitative  and  Quantitative  Analysis  of  Medieal  Ultrasound 


by 


Kevin  W.  Ayer 

Major,  United  States  Air  Force 


B.S.,  Chemistry,  University  of  Texas  at  El  Paso,  El  Paso,  Texas,  1982 
B.S.,  Electrical  Engineering,  Auburn  University,  Auburn,  Alabama,  1986 
M.S.,  Electrical  Engineering/Electro-Optics,  Air  Force  Institute  of  Technology, 
Wright-Patterson  Air  Force  Base,  Ohio,  1989 


A  dissertation  submitted  in  partial  fulfillment  of  the 
Requirement  for  the  degree  of  Doctor  of  Philosophy 
To  the  Chester  F.  Carlson  Center  for  Imaging  Science 
Rochester  Institute  of  Technology 


20001013  056 


CHESTER  F.  CARLSON  CENTER  FOR  IMAGING  SCIENCE 
ROCHESTER  INSTITUTE  OF  TECHNOLOGY 
ROCHESTER,  NEW  YORK 

CERTIFICATE  OF  APPROVAL 


Ph.D.  DEGREE  DISSERTATION 


The  Ph.D.  Degree  Dissertation  of  Major  Kevin  W.  Ayer 
has  been  examined  and  approved  by  the 
dissertation  committee  as  satisfactory  for  the 
dissertation  required  for  the 
Ph.D.  degree  in  Imaging  Science 


N.  A.  H.  K.  Rao,  Ph.D.,  Dissertation  Advisor 


Lb. 


Theodore  W.  Wilcox,  Ph.D. 


Date 


DISSERTATION  RELEASE  PERMISSION 
ROCHESTER  INSTITUTE  OF  TECHNOLOGY 
CHESTER  F.  CARLSON  CENTER  FOR  IMAGING  SCIENCE 


Title  of  Dissertation: 

Nonlinear  Ultrasoimd  Image  Modeling:  Development  of  a  complete  End-to-End  Model  for 
Qualitative  and  Quantitative  Analysis  of  Medical  Ultrasound 


I,  Kevin  W.  Ayer,  hereby  grant  permission  to  Wallace  Memorial  Library  of  R.LT.  to  reproduce 
my  thesis  in  whole  or  in  part.  Any  reproduction  will  not  be  for  commercial  use  or  profit. 


u 


Nonlinear  Ultrasound  Image  Modeling:  Development  of  a  complete  End-to-End  Model  for 
Qualitative  and  Quantitative  Analysis  of  Medical  Ultrasound 

by 

Kevin  W.  Ayer 

Submitted  to  the  Chester  F.  Carlson  Center  for  Imaging  Science 
in  partial  fulfillment  of  the  requirement  for  the 
Doctor  of  Philosophy  Degree 
at  the  Rochester  Institute  of  Technology 

Abstract 

Finite  amplitude  sound  propagation  imdergoes  nonlinear  distortion  due  to  continuous  path 
interaction  with  the  propagation  medium.  This  distortion  tends  to  defocus  the  beam  causing 
significant  lateral  and  contrast  resolution  degradation.  Fxmdamental  understanding  of  this 
interaction  requires  development  of  computational  models  that  accurately  predict  the  nonlinear 
interaction  —  development  of  media-bome  harmonics  —  as  well  as  produce  an  ultrasound  image  — 
introduction  of  transducer  effects,  interface  transitions,  and  innovative  image  processing  to  extract 
harmonics.  Most  computational  models  of  ultrasound  propagation  assume  axial  symmetry  for 
computational  expediency.  Two  notable  exceptions  are  the  K-Z-K  and  NLP  models.  A  new  end- 
to-end  model,  NUPROP,  is  introduced  that  also  incorporates  non-axially  symmetric  geometries 
and  simplified  transducer  responses  to  accurately  predict  ultrasound  RF  signals  for  image 
reconstruction.  Nonlinearities  are  modeled  using  either  the  Fubini  solution  or  Burgers’  Equation 
coupled  with  angular  spectrum  propagation  or  Lommel  formulation,  appropriately  masked  by  the 
transducer  frequency  response.  Comparative  analyses  are  performed  on  NUPROP  results  with 
high  correlation  with  the  literature.  Parameter  sensitivity  analyses  are  performed  to  determine 
harmonic  signal  characteristics  as  a  function  of  propagation  distance.  A-line  and  B-scan  images 
are  produced. 
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size  of  the  computational  space  allocated  for  signal  and  image  generation.  Environmental 
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Figure  6-8.  NUPROP  Tissue  Generation  Module.  Settings  used  to  generate  simulated  ultrasound 
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45  degree  angle  (black).  Object  3  is  the  inner  ellipse  (gray) . 64 
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6-7  are  redisplayed.  The  histogram  of  the  fundamental  image  differences,  second  harmomc 
image  differences,  and  third  harmonic  image  differences  are  presented.  The  third  harmonic 
difference  image  histogram  shows  appreciable  distribution,  whereas  the  fimdamental  and 
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Chapter  1*  introduction 

Medical  ultrasound  has  been  a  mainstay  in  noninvasive  imaging  of  the  human  body. 
Ultrasound  has  proven  to  be  a  relatively  risk  free,  accurate  and  versatile  technique  for  clinical 
diagnostics.  Ultrasound  allows  for  real-time  images  of  moving  structures  and  is  relatively 
inexpensive,  both  features  that  make  ultrasound  a  desirable  imaging  modality.  Even  though 
ultrasound  provides  such  benefits  for  medical  imaging,  two  reasons  keep  ultrasound  from  being 
more  widely  used  for  quantitative  medical  diagnostics;  first,  the  interaction  of  sound  with 
biological  tissue  is  a  very  complex  process,  with  many  parameters  that  are  either  assumed 
constant  or  measured  in  a  coupled  way.  A  coupled  measurement  implies  that  no  single  parameter 
can  be  measured  without  disturbing  or  affecting  the  measure  of  another  parameter.  For  instance, 
the  measure  of  absorption  depends  on  acoustic  impedance  or  scattering,  in  turn,  all  three  are 
dependent  on  temperature  and  temporal  frequency  of  the  input  signal.  The  second  reason  is  our 
lack  of  knowledge  of  fimdamental  interactive  processes.  One  such  interaction  results  in  the 
generation  of  harmonics  and  sub-harmonics  of  the  fimdamental  temporal  frequency  of  the  input 
signal.  This  phenomenon  is  known  as  nonlinearity. 

Harmonic  energy  generation,  or  nonlinear  wave  propagation,  is  frequency  dependent,  signal 
amplitude  dependent,  and  is  highly  susceptible  to  media  attenuation,  either  through  absorption  or 
scattering.  It  is  this  susceptibility  and  subsequent  lack  of  experimentally  derived  evidence  that 
gave  rise  to  the  false  assumption  that  nonlinear  effects  were  not  a  significant  contributor  to 
ultrasoimd  diagnostic  imaging.  Recent  advances  in  digital  architectures,  increased  dynamic  range, 
and  improved  signal  processing  capabilities  make  it  possible  to  exploit  harmonic  or  nonlinearly 
generated  energy  for  tissue  diagnostics.  This  thesis  provides  both  a  theoretical  underpinning  and 
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analytical  tool  to  understanding  how  nonlinear  ultrasound  waves  are  formed,  how  they  propagate, 
and,  more  importantly,  how  they  may  interact  with  biological  tissue  to  produce  significantly 

improved  ultrasound  images. 

Experimental  verification  and  validation  of  the  analytical  tool  (called  JVonlinear  Wtrasound 
PiJOPagation:  NUPROP)  is  presented.  Results  fi-om  exercising  NUPROP  illustrate  how  this 
analytical  tool  can  be  used  to  predict  fundamental  harmonic  wave  generation  properties  as  well  as 
investigate  the  influence  various  systematic,  environmental,  and  object  dependent  parameters  have 
on  image  formation. 
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Chapter  2.  Literature  Review 

This  chapter  discusses  the  current  literature  associated  with  nonlinear-linear  ultrasound  wave 
propagation,  namely  the  generation  of  nonlinear  longitudinal  waves  using  either  Burgers  or 
Fubini’s  approaches,  another  nonlinear  technique  using  the  K-Z-K  equation,  linear  wave 
propagation  using  angular  spectrum  or  the  Lommel  formulation,  and  image  formation  through 
modeling  of  the  imaging  system  as  a  linear,  separable  function  convolved  with  a  continuum  of 
point  scatterers  describing  the  medium  property  of  interest.  The  sections  are  divided  into 
Nonlinear  Wave  Propagation,  Linear  Wave  Propagation,  Other  Approaches,  and  Image 
Formation. 

Nonlinear  Wave  Propagation 

In  his  1997  paper,  Christopher  discusses  the  use  of  computational  predictions  based  on  his 
earlier  work,  and  their  application  to  experimental  validation,  whereby  second  harmonics  are 
produced  by  finite  amplitude  distortion.  Christopher  employed  a  fi-equency  domain  solution  to 
Burgers’  equation  to  derive  the  nonlinear  propagation  behavior  of  the  medium.  Burgers’  solution 
is  mathematically  developed  in  Appendix  I  and  described  in  Chapter  3  as  it  applies  to  this 
research. 

Christopher  described  his  experiment  wherein  he  used  a  concentric  ring  transducer  introducing 
a  sonic  waveform.  The  central  portion  of  the  concentric  design  used  a  circular  disk  operating  at  a 
5.0MHz  center  fi-equency.  The  concentric  ring  surrounding  the  disk  operated  at  2.5  MHz  and 
introduced  said  fi-equency  tone  at  1.6  microsecond  bursts.  Christopher  simulated  a  two-way  2.5 
MHz  input  burst  signal  and  a  second  harmonic  (5.0  MHz)  propagated  4.4  cm,  with  corresponding 
integrated  profiles.  Christopher  found  that  the  nonlinear  parameter,  ft,  was  not  significant  to  the 


4 


results  of  his  simulations.  Christopher  then  experimentally  obtained  comparable  results  for 
contrast  resolution  for  a  scatterer  free  or  cyst  region  of  tissue.  His  results  suggest  that  a  single 
pulse  methodology  is  suflBcient  to  produce  viable,  second  harmonic  generation  images. 

Krishnan,  Hamilton,  and  O’Donnell  [1997]  reported  that  harmonics  generated  by  nonlinear 
wave  propagation  in  the  tissue,  the  same  harmonics  that  reduce  contrast  in  perfiised  agent 
experiments,  could  be  modeled  by  a  Taylor  series  expansion  of  the  equation  of  state  for  plane 
waves  in  liquid  media  under  adiabatic  conditions.  Their  approach  directly  applied  Christopher  and 
Parker  [1991]  mathematical  concepts  and  laboratory  experimentation  to  modeling  contrast  agent 
image  enhancement.  Their  approach  to  describing  the  nonlinear  effect  employed  the  nonlinear 
ratio,  B/A.  ^  and  5  are  the  coefficients  associated  with  the  first  two  terms  of  a  Taylor  series 


expansion  of  the  state  equation.  The  equation  of  state  is  given  as 


/ P-Po  ,  B  p-po 
p-p„=A  — —  +-  — — 

\  Po  j  Po 


where  p  and  po  are  pressure  and  static  pressure,  respectively  and  p  and  po  are  density  and  static 


density,  respectively. 

Krishnan,  et.  al,  [1997]  developed  a  method  of  studying  harmonic  cancellation,  whereby  the 
second  harmonic  was  obtained  through  the  use  of  a  frequency  solution  to  Burgers’  equation 
(Haran  and  Cook  [1983];  Trivett  and  Buren  [1984]).  Krishnan,  et.  al,  [1997]  used  a  incremental 
step  approach  whereby  the  wavefront  was  propagated  a  distance,  Zi,  by  performing  a  linear 
propagation  foUowed  by  a  nonlinear  propagation.  The  linear  propagation  accounted  for  each 
frequency  component  while  the  nonlinear  propagation  accounted  for  the  generation  of  harmonics 
and  sub-harmonics.  Linear  propagation  was  accomplished  by  decomposing  the  signal  into  an 
angular  spectrum  of  plane  waves  (described  later  in  the  section  on  angular  spectrum).  Nonlinear 
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propagation  was  accomplished,  as  previously  mentioned,  using  a  frequency  domain  solution  to 
Burgers’  equation.  The  solution  takes  on  the  following  form 

.  ^TlftsZ 

(z  +  Az, /)  =  (z  +  Az, 0  +  y  2^2 

( n-\  ^  ^ 

2  ku'k  (z  +  Az,  i)u„_^  (z  +  Az,  0  +  2  mUz  +  O"*-™  (z  +  Az,  /)*  , 

n  =  1,2,. ..,A^. 

where  «„(z  + Az,/)  is  the  «'*  temporal  frequency  component  of  the  normal  velocity  waveform  at 
the  /'*  azimuthal  point  in  space  at  a  depth,  z  +  Az,  and  m„(z  +  Az,/)  is  the  same  temporal 
frequency  component  of  the  waveform  prior  to  nonlinear  propagation.  The  variable,  is  a 

function  of  the  nonlinear  ratio,  B/A,  as  =  1  +  — — .  The  variable  f  is  the  center  frequency  of  the 

2A 

propagating  waveform  and  c  is  soimd  speed  in  the  medium.  Krishnan,  et.  al,  did  not  consider 
attenuation  in  their  development. 

Another  technique  for  describing  nonlinear  wave  behavior  was  developed  by  Fubini  [1935], 
where  he  derived  an  explicit  solution  to  the  equation  for  motion  of  plane  soimd  waves  in  a  gas. 
Derivation  of  the  Fubini  solution  is  given  in  Appendix  II  and  described,  in  greater  detail,  in 
Chapter  3  relative  to  its  application  in  this  research.  Important  features  of  the  Fubini  solution  are 

1.  A  slope  term  in  the  denominator  that  reaches  negative  infinity  at  a  finite  propagation 
distance  known  as  the  shock  distance.  Computation  beyond  this  point  produces  invalid 
results 

2.  The  fimdamental  component  signal  amplitude  of  the  wave  will  decrease  as  the  wave 
propagates,  harmonic  component  signal  amplitudes  increase 

3.  As  initial  wave  amplitude  increases,  shock  distance  will  decrease 
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4.  As  frequency  increases,  the  shock  distance  will  decrease 

5.  The  ratio,  u/uo,  can  be  viewed  as  appropriately  weighted  sum  of  frequency  dependent, 
sinusoidal  basis  functions 

Linear  Wave  Propagation 

Lewin  [1989]  provided  a  more  detailed  description  of  the  angular  spectrum  technique.  His 
aim  was  to  provide  an  analytical  means  of  predicting  transducer  performance.  His  assessment  of 
the  advantages  of  using  the  angular  spectrum  technique  included  high  spatial  resolution, 
computational  efficiency  using  the  FFT,  and  the  method’s  applicability  to  both  forward  and  back 
propagation.  Lewin’s  derivation  of  the  angular  spectrum  method  is  repeated  in  the  following: 
Note:  this  derivation  can  be  found  in  several  references  Goodman  [1968];  Stepanishen  and 
Benjamin,  [1982];  Powers  [1976]: 

A  complex  pressure  wave  field  is  generated  by  a  set  of  monochromatic  sources.  A  wave  field 
is  represented  as  M(x,y,zJ  in  a  plane  paraUel  to  the  X-Y plane  at  a  position  as  depicted  in 

Figure  2-1. 


Figure  2-1.  Graphical  representation  of  a  complex  pressure  field  propagating  fi-om  one  parallel  plane  to  another  using  the 
angular  spectrum  technique. 
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The  two-dimensional  Fourier  transform  of  this  complex  pressure  wave  field  is  represented  as 
U(^Tj,Zo).  Each  complex  field  fi-equency  is  modified  by  an  appropriate  phase  fector  as  it 
propagates  fi'om  its  original  position  (zo)  to  a  new  position  (zj).  This  accounts  for  a  plane  wave 
propagating  fi-om  Zo  to  z;.  The  phase  factor  can  be  represented  as 


G(^,rj,z„,z^)  =  e 


(2-3) 


Consequently,  the  field  at  Z/  can  be  represented  as 
tt(jc,  7,  z, )  =  7,  z„ )  • 

where  A  is  the  temporal  wavelength  of  the  propagating  wave  field,  4^  and  7  are  the  angular  spatial 
fi-equencies  associated  with  the  x  and  y  coordinates,  respectively,  and  U(^tj,Zo)  is  the  Fourier 
transform  of  the  complex  field  at  z©.  Essentially,  the  propagation  of  a  complex  field  at  Zo  to  a 
position,  zj,  can  be  modeled  as  a  space-invariant  fiJ^ter  in  the  frequency  domain.  Notably,  if  /I  ^  is 
greater  than  ^  +  rf,  then  the  phase  term  is  purely  imaginary  giving  rise  to  harmonic  wave 
propagation.  Whereas,  if  the  converse  is  true,  then  the  phase  factor  is  real,  yielding  strong 
attenuation  in  the  positive  z-direction. 

Lewin  handled  evanescent  waves  in  the  back-propagation  direction  by  nulling  them  out  prior 
to  back-propagation  calculations,  thus  the  evanescent  waves  had  no  contribution  to  the  return  or 
reflected  signal.  Lewin  also  assumed  that  the  medium  acted  on  the  propagating  wave  the  same 
irrespective  of  propagation  direction.  Levrin’s  application  of  the  angular  spectrum  technique  was 
strictly  linear. 

Daly  and  Rao  [1998]  proposed  a  closed  form  fi-equency  domain  formalism  whereby  spatially 
integrated  diffraction  corrections  were  employed.  The  Lommel  diffraction  formulations  were 
used  as  a  “coupled-pair”  or  approximate  Fourier  transform  pairs,  amenable  to  closed-form  spatial 
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integration,  deriving  an  exact  diffraction  solution  for  a  circular  piston  transducer  geometry.  As 
Daly  and  Rao  point  out,  diffraction  corrections  are  necessary  when  trying  to  characterize 
biological  tissue  using  ultrasound,  especially  in  the  near  field,  sometimes  inappropriately  referred 
to  as  the  Fresnel  region.  Lommel  fimctions  are  defined  as  follows: 


/  \n+2s 

'  U  ' 


5=0 


Vvj 


f;(“.v)=£(-ir(«+25) 


.  n+25 


5=0 


\u) 


(2-5a) 

(2-5b) 


Since  the  Lommel  fimctions  are  infinite  summations,  they  can  only  be  approximated  numerically. 
In  the  Fresnel  approximation,  and  imposing  circular  symmetry,  the  velocity  potential. 


/f,  {p,  z,(o) ,  can  be  approximated  as 


^  \  -ik(z+^)  'j*  ik^ 

H^{p,z,co)  =  -e 

Z 


(2-6) 


where  p^  =  +  y]  is  the  off-axis  distance  at  z  =  0,  plane,  p  =  is  the  off-axis  at  the 


z  plane,  a  is  the  aperture  diameter,  and  the  hat  indicates  estimation.  Invoking  the  singularity 
function  with  a  change  of  upper  limit  integration  (equal  to  infinity),  then  (2-6)  can  be  interpreted 
as  the  Hankel  transform  of  the  product  of  the  singularity  function  and  a  quadratic  phase  term, 
rewriting  (2-6)  as 


1  -ik{z+^) 

H^{p,z,G))  =  -e 
Z 


—  J 
P 


\  z 


I 


(2-7) 


where  the  convolution  is  now  with  respect  to 


with  the  closed  form  solution  given  as 


1 

Hi{p,z,(o)  =  je 


-/  kz+ — +— 
2u  2 


[U,{u,v)  +  iU,iu,v)] 


(2-8) 
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with  u^ka^  I zsndv-kap! z . 


Oberhettinger  [1961]  developed  exact  expression  for  the  inverse  Fourier  transform  of  Hi 


given  as 
forp  <  a 


hiip,z,t)  = 


0, 

c. 


c 

— arccosi 
n 


{c(f  -z^  +  p^  -a^ 
2p4ictf  -z^ 


0, 


ct  <z 
z  <  ct ,  R' 

,  R'<ct<R 

ct>R 


for p>  a 


hy{p,z,t)  =  \ 


—arccosi 
n 


{ctf  -z^  +  p^  -  a 

2p^icty-z^ 


0, 


ct<R' 

,  R'<cl<R 
ct>R 


(2-9) 


(2-10) 


where  R'=  ^Jz^  +(a- pf  and  R  =  4z^  +{a  +  pf  .  With  significant  manipulation  and 

invocation  of  relationships  developed  by  Wolf  [1953];  Daly  [1998],  the  spatially  integrated, 
estimation  of  the  Fourier  domain  representation  of  the  velocity  potential  is  given  as 

(2-11) 


Izjf  \\  /•"  '\ 


where  u  =  ka^jz ,  a  Fresnel  distance  normalized  wavenumber. 


Other  Approaches 

Nikoonahad  and  Iravani  [1989]  incorporated  a  complex,  fi-equency  dependent  compressibility 
into  the  wave  equation  to  predict  the  focal  plane  distribution  in  biological  media.  Their  approach 
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employs  a  complex,  frequency  dependent  wavenumber  derived  from  ultrasonic  propagation  given 


by 


W^u-pk 


rd^u 


dt^ 


0 


(2-12) 


where  it  is  a  complex  wavenumber  defined  by  O’Donnell  [1981]  as 


k  =  kj((D)  +  ik2(0)) 


A:,(o)  = 


1 


C(0)) 


k^ico)  = 


2a(<y) 

poacico) 


y 

) 


(2-13) 


where  c(c^  and  a(a^  are  phase  velocity  and  absorption  coefficient,  respectively.  The  subscript  on 
c  refers  to  the  phase  velocity  at  reference  frequency,  fit-  Nikoonahad  and  Iravani  state  that  the 
absorption  coefficient  must  exhibit  two  properties: 

1 .  Nearly  linear  dependence  on  frequency 

2.  A  square  law  dependence  on  frequency  in  liquids  where  relaxation  times  are  short. 
Nikoonahad  and  Iravani  reported  that  Seghal  and  Greenleaf  [1982]  used  these  two  properties 


and  formulated 


f  c(fi))^ 

1  J 

n 

[^(fi>,)  +  5(fi>,)J 

a{o})  = 


2fi!^(®)ln(y) 


c(fi))5(fi)) 
=  <ur 


B(ci))  =  yll  +  icory 


(2-14) 


For  these  relationships,  yis  the  specific  heat  ratio  and  r  is  the  relaxation  time.  By  substituting 
these  relationships  into  the  previous  equation,  the  following  relationships  for  ki  and  are  derived 
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(2-15) 


A:,(«)  = 


PCr 


A{G),)  +  Bi(0,) 
A(co)  +  B((i}) 


(4ln/y 


kA(o)  = 


4>4(ftj)  In  X  r  MO>r  )  +  ) 

pclBio))  1,  A{a))  +  B{a}) 


Both  ki  and  itj  accovint  for  dissipative  (absorption)  and  dispersive  (frequency  dependent 
attenuation)  processes  within  biological  media.  It  is  this  complex  wavenumber  that  is  introduced 
into  a  three-dimensional  Fourier  transform,  where  the  angular  dimensions  are  x  and  y,  and  the 
third  dimension  is  time,  t.  Nikoonahad  and  Iravani  then  employed  the  traditional  angular 
spectrum  technique  to  linearly  propagate  the  transverse  wavefront,  using  the  complex,  frequency 
dependent  wavenumber  as  a  modifier  fi)r  the  propagation  term.  Unfortunately,  this  approach  is 
strictly  linear  and  does  not  account  for  harmonic  component  generation. 

Averkiou,  Roundhill,  and  Powers  [1997]  proposed  a  different  mathematical  approach 
employing  the  Khokhlov-Zabolotskaya-Kuznetsov  (KZK)  equation  developed  by  Zabolotskaya 
and  Khokhlov  [1969];  Kuznetsov  [1970].  The  KZK  relationship  accounts  for  diffraction, 
absorption,  and  nonlinearity.  Averkiou,  et.  al.,  concluded  that,  since  tissue  properties  were 
readily  available  in  literature,  including  coefficients  for  nonlinearity,  this  tissue  property  could  be 
modeled  using  the  KZK  equation.  The  KZK  model  is  given  as 

n  C  S  D  B  (2"16) 

- -2-\7^  n -i-  —  H - — - — 

dzdt  2  2cl  dt^  2p„cl  dt^ 


where  p  is  the  sound  pressure,  z  is  the  coordinate  along  the  propagation  direction. 


=  r  is  the  angular  radial  coordinate  (assuming  axial  symmetric  sound  beam), 

"  dr^  rdr 

t  =  t'-—  the  retardation  time,  and  is  the  sound  speed.  The  first  term  on  the  right-hand  side 
Co 
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of  this  relationship  accounts  for  diffraction;  the  second  term  accounts  for  thermoviscous 
attenuation  or  dissipation  where  5  is  the  diffiasivity  of  sound  LighthiU  [1980];  the  third  term 

accounts  for  nonlinearity,  where  p  is  given  as  1  +  ^  Beyer  [I960].  One  limitation  of  this 

approach  is  the  requirement  for  a  directive  sound  source  where  ka»\,k  bemg  the  wavenumber 
associated  with  the  source  and  “a,”  the  radius  of  the  source.  Another  limitation  is  the  paraxial 
approximation,  thus  limiting  the  vaUdity  of  solutions  to  regions  around  the  axial  direction. 
Although  the  KZK  equation  is  not  universal,  it  does  provide  good  agreement  with  experiments 
adhering  to  these  restrictions  Averkiou  and  Hamilton  [1997];  Baker  and  Humphrey  [1990], 
Baker,  Berg,  and  Tj<|)tta  [1997]. 

Image  Modeling 

Image  formation  is  generally  divided  into  two  main  tasks;  modeling  the  imaging  system  and 
modeling  the  tissue.  Most  simulations  assume  a  constant  velocity  throughout  the  propagation 
medium,  zero  attenuation,  zero  noise  level,  infinite  receiver  bandwidth,  and  linear  gain 
characteristics  of  the  receiver.  An  exceUent  example  of  this  type  of  simulation  was  developed  by 
Bamber  and  Dickinson  [1980].  Bamber  assumed  that  the  system  forming  the  image  had  a  linear, 
separable,  space  invariant  point  spread  fimction,  with  the  tissue  or  other  object  represented  by  a 
continuous  distributed  set  of  point  scatterers.  Bamber’s  assumptions  then  yielded  an  image, 
I(x,y),  using  the  convolution  (denoted  ®)of  the  system  point  spread  function,  H(x,y)  and  the  set 
of  point  scatterers,  T(x,y) 

I{x,  y)  =  H{x,  y)  <8>  r  (x,  y)  (2-17) 

Pulse  echo  imaging  is  a  tomographic  technique  that  differs  fi-om  conventional  imaging  in  a  number 
of  respects:  signal  beam  is  coherent  and  the  detector  is  phase  sensitive.  This  will  significantly 
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impsct  the  iniage  reconstruction  since  return  signuls  from  diflfering  locations  will  interfere,  giving 
rise  to  classical  speckle  effects.  Bamber’s  model  neglected  multiple  scattering.  Additionally, 
Bamber’s  model  assumed  non-axial  symmetry  such  that  the  point  spread  function  consists  of  two, 
separable  convolution  filters:  one  for  the  longitudinal  pulse  shape  and  the  other  for  the  trans-axial 
beam  profile  such  that  H(x,y)  firom  equation  2-17  becomes  H 2  (x,l)  and  //,  (l,y)  yielding 

Iix,y)  =  H2ix,l)  <S)  H,  (l,y)  0  r(x,y)  (2-18) 

Convolution  with  H2(x,\)  equated  to  detection  if  x  was  the  longitudinal  direction  while 
convolution  with  Hi(l^)  equated  to  smoothing  or  diffraction  if y  was  the  trans-axial  direction. 

As  previously  mentioned,  Bamber’s  model  depended  on  describing  the  medium  as  a 
continuous  distribution  of  point  scatterers,  or  the  acoustic  impulse  response  of  the  tissue.  One 
popular  tissue  model  used  by  Bamber  was  the  inhomogeneous  continuum  models  of  Chivers 
[1977],  Gore  and  Leeman  [1977],  and  Nicholas  [1977]  in  which  the  density,  p,  and 
compressibility,  of  the  tissue  were  modeled  as  continuous  and  fluctuated  about  their  mean 
value.  Fields  and  Dunn  [1973],  however,  neglected  density  fluctuations  opting  for  compressibility 
fluctuations  only.  Bamber  suggested  that,  if  the  medium  was  considered  comprised  of  a 
continuum  of  elemental  scatterers  within  an  inhomogeneous  medium,  and  the  density  term  being 
neglected,  then  the  source  strength  of  each  elemental  scatterer  was  proportional  to  the  second 
derivative  of  the  compressibility  at  each  elemental  position  in  the  medium.  This  is  an  important 
concept  because,  if  the  spatial  distribution  of  the  tissue  compressibility  is  known,  then  the  required 
tissue  impulse  response  can  be  derived  such  that 

Tix,y,z)  =  V^P{x,y,z)  (2-19) 
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Bamber  further  stipulated  that  P(x,y)  had  to  be  purely  real,  constraining  the  real  spectrum  to 
be  an  even  function  of  frequency  and  the  imaginaiy  spectrum  to  be  an  odd  function  of  frequency. 
After  obtaining  the  compressibility  matrix  values,  Bamber  then  formed  images  by  sunply 
multiplying  the  Fourier  domain  representation  of  the  compressibility  matrix  by  s^. 

Chapter  2:  Summary 

The  literature  is  replete  with  examples  of  ultrasound  modeling  techniques  that  capture  aspects 
of  the  overall  process.  For  nonlinear  wave  propagation,  two  models  are  promment; 
Christopher’s  and  Parker’s  NLP  and  the  Khokhlov-Zabolotskaya-Kuznetsov  (KZK) 
implementation.  Christopher’s  NLP  is  based  on  solutions  to  Burgers’  equation  and  is 
implemented  using  a  C-based  code.  Averkiuo’s  Implementation  of  the  K-Z-K  includes 
diffiaction,  attenuation,  and  nonlinearities  in  a  single,  third  order  partial  differential  equation. 

For  linear  propagation,  two  methods  are  presented:  angular  spectrum  and  Lommel 
formulation.  Angular  spectrum  implements  a  planar  geometry  derived  mathematical  convention 
using  a  propagation  factor  that  accounts  for  complex  field  changes  over  small  incremental 
distances.  Since  angular  spectrum  is  developed  using  planar  calculus,  focused  beam  geometries 
are  not  easily  accommodated. 

Angular  spectrum  is  applicable  throughout  the  propagation  path,  from  just  in  front  of  the 
transducer  to  as  far  as  the  user  wished  to  compute  the  field.  Additionally,  angular  spectrum  is 
suitable  for  computing  field  propagation  from  non-radially  symmetric  geometries. 

Lommel  formulation  implements  planar  as  well  as  focused  beam  transducer  geometries. 
However,  Lommel  is  appUcable  for  circularly  symmetric  geometries  only.  Additionally,  Lommel 
is  strictly  a  linear  propagation  process.  Therefore,  any  process  that  introduces  nonlinearities,  i.e. 
generation  of  harmonics,  will  necessarily  be  inaccurately  predicted  by  Lommel. 
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Ultrasound  image  generation  models  depend  strongly  on  the  tissue  characterization  or  tissue 
construction  model  used.  One  model  has  found  considerable  favor  in  literature:  Bamber  and 
Dickinson’s  continuous  array  of  point  scatterers.  If  the  tissue  characteristic  of  interest  can  be 
modeled  as  an  array  of  continuous  point  scatterers,  then  Bamber’ s  techmque  is  reasonably 
accurate  in  predicting  linear  ultrasound  effects.  Bamber’s  technique  involves  the  convolution  of 
the  tissue  array  with  the  point  spread  function  of  the  ultrasound  pulse.  Bamber  accomplishes  the 
convolution  by  taking  the  process  to  the  Fourier  domain  and  multiplying  the  Fourier  transforms  of 
the  tissue  array  and  the  point  spread  fimction,  then  inverse  Fourier  transforming  the  product. 

Though  models  exist  for  the  various  aspects  of  the  ultrasound  process,  no  single  model 
encompasses  all  aspects  of  ultrasound  propagation  and,  ultimately,  ultrasound  imaging.  In  this 
work,  a  new  model.  Nonlinear  Ultrasound  Propagation  (NUPROP),  is  introduced  that  captures 
aspects  of  nonlinear  wave  generation  and  propagation,  linear  wave  generation  and  propagation, 
and  image  generation. 


Chapter  3*  NUPROP  Module  Development 

This  chapter  discusses  the  main  modules  comprising  NUPROP:  1)  phantom  generation,  2) 
input  signal  waveform  type  (sinusoid  or  chirp),  3)  implementation  of  wave  propagation  - 
Nonlinear:  Burgers’  Equation  and  Fubini  Solution,  and  Linear:  Angular  Spectrum  and  Lommel,  4) 
image  generation  using  the  Hilbert  transform,  and  5)  post  processing  such  as  statistical  analyses. 

Phantom  Generation:  Continuous  Point  Scatterer  Array 

The  validity  of  simulating  ultrasound  images  relies  on  the  ability  to  describe  the  object  or 
tissue  as  a  continuous  array  of  point  scatterers:  a  model  of  the  acoustic  impulse  response  of  the 
tissue.  This  view  of  modeling  tissue  characteristics  is  supported  by  the  theory  of  scattering  from 
an  inhomogeneous  medium.  Phantoms  can  be  represented  mathematically  as 

Ph{x,  y,  z)  =  A{x,  y,z)- A  ^  ^ 

where  A(x,y,z)  is  the  amplitude  of  the  acoustic  impedance  at  a  point  in  space,  A  is  the  mean  value 
of  the  array,  vAih.  A(x,y,z)  computed  using 

A{x,y,z)  =  ^a„RANDOMU(seed,x„,y„,z„)  (3-2) 

n 

given  cca  is  the  object’s  acoustic  impedance  value  at  position  x„,  y„,  z„,  RANDOMU  is  a 
uniformly  distributed  random  number  generating  function,  and  seed  is  a  temporally  changing  value 
to  initiate  the  generation  of  a  random  number.  If  seed  is  not  specified,  then  the  algorithm  will  use 
the  computer’s  clock  to  generate  a  seed  value.  The  mean  value  is  subtracted  from  the  array 
composite  to  eliminate  any  artificial  discontinuities  introduced  by  the  phantom  generation 
algorithm. 
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Waveform  Type 


Single  Frequency  Sinusoid 

Since  a  single  frequency  is  equivalent  to  a  continuous  vfawe  input,  it  has  little  to  no  use  in 
imaging,  since  the  return  signals  cannot  be  distinguished  from  one  point  scattering  position  to 
another.  However,  the  single  frequency  sinusoidal  input  is  useful  in  comparing  NUPROP’s 
nonlinear  wave  propagation  performance  against  existing  data  in  literature  (Ryan  s  Fubini 
solution,  Christopher’s  trans-axial  beam  profiles,  and  Averkiuo’s  along-track  or  longitudinal  beam 
amplitude  profile). 

Narrowband  Frequency  Sinusoid 

The  narrowband  or  pulsed  sinusoidal  input  is  typically  used  ftir  image  generation  (B-scan). 
The  narrowband  sinusoidal  input  requires  care  in  constructing  the  pulse  or  temporal  apodization 
function.  Whether  the  pulse  apodization  is  frequency  dependent  or  not,  the  temporal  extent  must 
be  such  as  to  produce  minimal  frequency  domain  overlap  of  the  bands  about  the  harmonics. 
Figure  3-1  illustrates  this  condition.  Note  that  the  frequency  dependent  apodization  curves  have 
greater  overlap  than  the  frequency  independent  curves. 


Figure  3-1.  Frequency  dependent  and  independent  apodi2ation  schemes  for  ensuring  minimum  overlap  of  frequency  signals: 
fundamental,  second,  and  subsequent  harmonics.  The  signal  is  presented  in  dB  scale.  The  fundamental  frequency  is  arbitrarily 
set  to  2.5  MHz,  with  second  and  third  harmonics  shown  at  5.0  MHz  and  7.5  MHz,  respectively.  The  apodization  functions  are 
Gaussian  with  exponent  2. 
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Linear  Frequency  Modulated  Sinusoid  (Apodized  and  Non-Apodized) 

Another  multiple  frequency  signal  input  is  the  linear  frequency  modulated  sinusoid  more 

commonly  called  the  “chirp.”  Arguably  one  of  the  most  useful  linear  filtering  techniques  is  the  use 
of  the  quadratic  phase  filter  or  chirp.  The  primary  advantage  of  using  a  chirp  signal  over  a  single 
frequency  sinusoid,  or  apodized  sinusoid,  is  the  ability  to  have  a  very  wide  bandwidth,  infinite  in 
fact,  while  maintaining  reasonable  signal  gain  (or  signal  to  noise  ratio).  Figure  3-2  shows  both  the 
apodized  and  non-apodized  chirp  signal  for  start  frequency  of  1.5  MHz,  sweep  rate  of  8.5x10'®. 
The  apodization  used  for  the  chirp  ensures  minimum  overlap  of  frequency  harmonic  bands. 


Real  Imt^nary 


Figure  3-2.  Examples  of  apodized  and  non-apodized  linear  frequency  modulated  sinusoids.  Initial  frequency:  1.5  MHz,  sweep 
rate  of  1.5x10^^  and  apodization  to  ensure  minimum  harmonic  signaJ  overlap.  Note  that  the  chirps  represented  here  consist  of 
the  fundamental  and  1 0  harmonics  of  the  initial  frequency. 

The  non-apodized  chirp  input  signal  has  the  same  diagnostic  advantages  as  the  single 
frequency  sinusoid,  but  suffers  the  same  non-imaging  capability.  The  apodized  chirp  possesses 
the  imaging  capabihty,  but  suffers  strict  apodization  requirements  to  ensure  minimum  overlap. 
These  restrictions  may  prove  to  be  too  stringent  so  as  to  negate  the  apodized  chirp’s  inherent 
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processing  advantage:  improved  signal  to  noise  ratio  without  necessary  bandwidth  reduction.  The 
development  of  NUPROP’s  implementation  of  the  chirp  is  given  in  Appendix  III. 


Wave  Propagation:  Nonlinear 

Burgers’  Equation  Model 

The  Burgers’  equation  models  the  effect  of  wave  shearing  in  viscous  and  non-viscous  media. 
The  Burgers’  relationship  can  be  represented  as  a  second  order,  non-linear  partial  differential 
equation  as  follows: 

du  du  „2  n 

dt  dz 

where  u  represents  a  complex  wave  field  and  v  is  the  viscosity  of  the  media.  Solutions  to  the 
Burgers’  equation  are  of  the  form  u  =  f(Qt-kz)  where  k  is  the  wavenumber,  z  is  the  propagation 
direction,  <o  is  the  wave  frequency,  and  t  is  time.  Blackstock,  as  referenced  in  Beyer  and  Letcher 
[1969]  developed  a  closed  form  solution  to  Burgers’  equation  given  as 

(3-4) 
u 


/F=l 

r  r 

aiUBITA)  1 
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where  a  is  the  effective  absorption  coefScient,  Co  is  the  wave  speed,  o)  is  the  wave  frequency, 
(1+B/2A)  is  the  nonlinear  parameter,  k  is  the  wavenumber,  /  is  the  shock  distance  (see  3-8),  or 
discontinuity  distance  as  described  in  the  Fubini  Solution  section  of  this  chapter,  and  In  is  the 

Bessel  function  of  order  n  of  the  imaginary  argument 

=  (3-5) 
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Fubini  Solution  Model 

Another  technique  for  describing  nonlinear  wave  behavior  was  developed  by  Fubini  [1935], 
where  he  derived  an  explicit  solution  to  the  equation  for  motion  of  plane  sound  waves  in  a  gas 

c] 

(I +  d^  I  day"'  da^ 


as 


00 
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(3-7) 
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A  discontinuity  will  occur  when  the  denominator  term  of  (3-6),  d^jda,  becomes  mfimtely 


negative.  This  “slope”  term  must  approach  negative  infinity  at  some  propagation  distance  due  to 
the  fact  that  large  values  of  u  will  propagate  faster  than  small  values  of  u,  thus  giving  rise  to  wave 
distortions.  This  phenomenon  is  commonly  referred  to  as  a  “shock  fi'ont.”  Fubini  computed  this 

distance  fi-om  the  origin  to  the  point  where  the  shock  forms  as 

^2  (3-8) 

/=— - - ^ - 

(,Bl2A  +  \)(ou„ 

Appendbc  II  describes  the  mathematical  development  of  this  closed  form  solution  to  Fubini’s 


equation. 

Up  to  this  point,  all  solutions  have  been  for  continuous  or  undamped,  single  fi-equency 
sinusoidal  input.  A  more  common  ultrasound  imaging  modality  is  the  pulsed  sinusoidal  signal, 
whereby  a  single  fi-equency  drive  signal  is  modified  using  an  envelope  or  temporal  windowing 
function.  Such  a  drive  signal  can  be  expressed  as 

/ (z,/)  =  giz,  t)  •  sin(6jr  -  kz)  (3-9) 

where  g(z,t)  describes  the  windowing  function.  A  typical  windowing  function  would  be  a 
Gaussian.  By  “windowing”  the  sinusoidal,  the  input  signal  is  no  longer  single  fiequency  but  rather 
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a  “band”  of  frequencies  centered  on  the  fimdamental.  If  more  than  one  frequency  exists  in  the 
drive  signal,  then  each  component  must  be  modified  by  an  appropriately  frequency  scaled 
Gaussian  to  account  for  this  frequency  “broadening.”  If  sufficient  separation  exists  between 
multiple  frequencies,  then  a  reasonable  approximation  to  the  overall  signal  would  be  a  simple 
modification  of  (3-7)  to  include  the  frequency  dependent  apodization  terms  such  that 

(3-10) 

where  a  is  an  engineering  scaling  fimction  to  ensure  sufficient  localization  of  energy  about  each 
frequency  term,  (t„  is  the  frequency  dependent  standard  deviation,  and  s  is  at-kz. 

Wave  Propagation:  Linear 

Angular  Spectrum  Model 

Linear  propagation  is  accomplished  in  NUPROP  one  of  two  ways:  Angular  Spectrum  or 
Lommel  Formulation.  This  section  discusses  NUPROP’s  implementation  of  Angular  Spectrum. 
As  previously  discussed,  angular  spectrum  is  a  technique  whereby  a  planar  wavefront  can  be 
translated  along  a  propagation  path,  usually  denoted  as  the  z-direction,  in  incremental  steps. 
NUPROP  assumes  the  pressure  amplitude  is  known  at  the  source  or  excitation  field  position 
designated  z  =  0.  Goodman  [1968,  ch.  3]  described  the  angular  spectrum  as  the  Fourier 
transform  of  the  pressure  field  with  respect  to  x  and  y,  given  values  of  the  pressure  field  in  the  x-y 
plane  at  a  distance,  z.  Defining  the  pressure  field  as  p(x,y,z),  then  its  angular  spectrum  is  given  as 

The  relationship  between  P(^ti>z)  and  P(^,tj,0)  comes  from  the  fact  that  p(x,y,z)  obeys  the 
Helmholtz  equation  (S7^  +  k^)p  =  0,  where  k  =  2^/ A  is  the  wavenumber,  resulting  in 


u  X-'  T  nz  \  . 

— =Z— T  r' 
Un  «=i  nz  y  I  J 


sin[n(5)]-e 
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P{^,  T],  z)  =  P{^, 

where  H(^,rj,z)  is  defined  as  the  space-invariant  propagation  filter 


(3-12) 


1 


-kzyli\p+n^)-i 


,forr+r  >-:;r 


(3-13) 


By  inspection,  if  +  ;7^  <  ^,  then  the  wave  propagates  with  an  amplitude  of  unity  in  the 


absence  of  attenuation.  If  +7^  >  ^,  then  the  wave  exponentially  decays  with  propagation. 

A, 


This  exponentially  decaying  wave  is  called  an  evanescent  wave.  Evanescent  waves  decay  within  a 
few  wavelengths  of  propagation  so  they  are  typically  ignored  for  simulation  purposes. 

Again,  by  inspection,  if  the  wavefi'ont  is  backpropagated  fi'om  an  array  of  known  data,  then 
the  evanescent  wave  will  grow  exponentially,  since,  by  convention,  the  sign  of  the  exponent  is 
positive.  Typically,  the  evanescent  wave  is  set  to  zero  prior  to  backpropagation  or  the  sign  of  the 
exponent  is  set  negative  to  account  for  evanescent  waves  decaying  in  actual  tissue.  NUPROP’s 
implementation  of  angular  spectrum  assumes  planar  apertures. 


Lommel  Formulation 

This  section  discusses  the  development  of  Lommel  formalism  relative  to  NUPROP’s 
implementation.  As  Daly  and  Rao  have  determined,  the  Lommel  formulation  is  a  closed-form 
estimate  of  the  Fourier  transform  of  the  arccos  diffraction  formulation  and  can  be  used  to 
determine  the  wavefi-ont  characteristics  at  a  propagation  distance  fi-om  an  extended  source.  The 
Lommel-based  formalism  allows  the  direct  calculation  of  diffraction  effects  in  the  fi-equency 
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domain,  in  the  absence  of  nonlinearities.  NUPROP  uses  Daly’s  IDL  code  to  compute  on-axis 
amplitudes  for  circularly  symmetric  apertures. 


Image  Generation 

The  wave  equations  used  for  generating  nonlinear  ultrasound  waves  (see  section  on  nonlinear 
wave  propagation)  represent  the  incident  wave.  Bamber  [1980]  computes  the  signal  ampUtude  at 
any  image  point  from  the  received  voltage  by  transforming  time  dependence  to  distance 
dependence,  deriving  the  equation  describing  the  image  formation  as 
I(x,y,z)=  \a(2x'-2x)b(y — y,z'-z)T(x',y',z')dh 


where  I(x,y,z)  is  the  image,  a  is  the  pulse  shape,  b  is  the  pulse-echo  beam  profile,  and  T  is  the 
tissue  impulse  response  given  as 


Adz\  p„  Po  j 


(3-15) 


where  pi(r)  is  the  spatial  density,  Po  is  the  static  density,  Pi(r)  is  the  spatial  compressibility,  and  Po 
is  the  static  compressibility.  This  relationship  holds  true  provided  the  wavelength  of  the  incident 
pulse  is  much  smaller  than  the  pulse  length  or  the  beam  width.  The  image  amplitude  distribution, 
then,  is  a  convolution  between  the  tissue  impulse  response  and  the  pulse-echo  point  spread 
function.  So,  NUPROP  computes  harmonic  component  contributions  at  each  penetration  depth 
using  nonlinear  processes  (longitudinal  pulse  profile),  coupling  those  processes  with  the  linear 
propagation  to  compute  the  angular  or  trans-axial  beam  profiles  at  those  penetration  depths,  then 
convolves  the  resulting  function  (a  three  dimensional  wave  function)  with  the  phantom  point 
scatter  array  at  each  position.  The  result  of  the  convolution  produces  the  RF  signal  for  a  given 
simulation.  B-Scan  images  are  then  generated  using  conventional  envelope  detection  via  the 
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Habert  transform.  Typically,  ultrasound  is  used  clinically  in  a  pulsed  mode,  a  short  duration 
temporal  “burst”  of  ultrasonic  energy.  The  pulse  then  travels  a  distance  into  the  medium,  reflects 
off  the  object  of  interest,  returning  to  the  same  transducer  that  produced  the  pulse.  The 
transducer  is  now  “switched”  to  receive  mode,  detecting  the  return  echoes  from  reflections  within 
the  medium.  The  time  taken  between  initial  pulse  and  return  signals  constitutes  the  collected  data. 
The  intensity  of  the  echoes  is  also  recorded,  and  is  used  when  rendering  the  scans. 

The  pulse  returns  are  converted  into  weak  voltage  signals  proportional  to  the  echo  intensity. 
The  converted  echo  signals  are  conditioned  through  a  series  of  electronic  circuits. 

1.  limiter  designed  to  prevent  excessively  large  transient  voltage  signals  to  damage  the 
receiver  electronic  circuitry, 

2.  a  logarithmic  amplifier  to  boost  the  weak  signal  amplitudes  (as  much  as  100  dB), 

3.  a  rectifier  to  convert  negative  half-cycles  of  the  echo  signals  into  positive  half-cycles, 

4.  a  demodulator  to  remove  the  carrier  frequency,  leaving  only  the  magmtude  envelope  of 
the  echo  signal. 

Envelope  detection  is  typically  performed  on  the  fundamental  frequency  earner.  The  demodulator 
performs  a  Hilbert  transform  on  the  data  to  produce  the  envelope  detection.  The  Hilbert 
Transform  provides  a  90  degree  phase  shift  to  its  input.  If  x(t)  is  the  real  input  signal,  and  the 
Hilbert  transformed  signal  is  h(t),  then  the  analytic  signal  is: 

z{t)  =  x(t)  +  ih(t)  (3*16) 

If  the  analytic  signal  possesses  a  single  frequency,  then  the  amplitude,  A(t),  of  the  waveform 
can  be  represented  as  the  square  root  of  the  squared  magnitude  such  that 

A(t)  =  yjz(t)  •  z(ty  =  [(x(/)  +  ihit)) •  (x(0 - =  ^Jx(ty  +h(tf 
where  z(t)*  represents  the  complex  conjugate  of  z(t). 
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Figure  3-3  graphically  illustrates  NUPROP’s  modular  design  and  processes  for  simulating 
nonlinear  ultrasound  waves  and  image  generation. 


Figure  3-3.  Modular  design  and  processes  associated  with  NUPROP:  Nonlinear  Ultrasound  Propagation  Modeling.  A) 
Transducer  Effects  and  Input  Waveform  Generation.  B)  Nonlinear  -  Linear  Coupled  Wave  Propagation.  C)  Wave  -  Object 

Interaction  and  Phantom  Generation.  D)  Linear  Propagation.  E)  Transducer  Effects  and  Image  Simulation. 


A I  introduces  transducer  effects  to  the  input  signal.  The  transducer  is  modeled  as  function 
that  introduces  beam  spatial  and  temporal  shaping.  The  spatial  beam  shaping  is  simple:  circular 
or  rectangular  geometry,  unit  amplitude  within  the  boimdaries  defined  by  the  geometry,  non- 
focused.  The  temporal  beam  shaping  is  simple  as  well;  Gaussian  with  exponent  varying  from  two 
to  sbc.  Gaussian  temporally  apodizing  functions  with  exponent  of  3  or  more  are  called  “Super- 
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Gaussian”  and  are  characterized  by  sharper  fall-offs  approaching  the  Rect  function  defined  as 


follows: 


1  for  <  0.5 

f(t)  =  RECT{t)  =  lim  =  ]  0.5  for  |/|  =  0.5 

0  otherwise 


(3-17) 


A2  depicts  the  pulse  ‘Vave-packet”  introduced  into  the  system  by  the  transducer.  Note  that 
the  along-track  profile  is  Gaussian  with  a  particular  carrier  fi-equency.  The  cross-track  profile  is  a 
“Super-Gaussian”  with  eiqionent  value  of  6.  B  shows  the  “wave-packet  propagating  through  the 
medium.  The  medium  is  assumed  to  be  homogeneous  throughout  the  propagation  path.  The  plot 
in  B  shows  the  along-track  profile  “sUce”  indicating  the  generation  of  harmonics  by  the  medium 
manifested  as  a  “sawtooth”  waveform  C  depicts  interaction  of  the  nonlinear  wave  vsdth  the 
object.  Note  that  the  temporal  signal  shows  significant  harmonic  contributions  to  the  waveform 
(sawtooth  nature  of  the  along-track  wave  profile).  Also,  diffraction  has  changed  the  cross-track 
profile  to  be  more  and  more  like  a  “sombrero”  function  for  circular  transducer  geometries.  The 
sombrero  function  is  defined  in  Easton  [1997]  as  a  radially  symmetric  function  proportional  to  the 


where  Jiinr)  is  the  Bessel  function  of  the  first  kind  of  order  unity.  The  Fourier  domain 
representation  of  the  spatio-temporal  wave  is  depicted  in  the  lower  portion  of  C.  Note  that  the 
fundamental  and  its  harmonics  are  present  in  the  signal,  as  indicated  by  signal  strength  at  the 
harmonic  frequencies.  After  the  nonlinear  wave  interacts  with  the  object,  reflection,  or  echo, 
signals  are  generated.  These  echo  signals  propagate  back  toward  the  sensing  transducer.  This 
propagation  path  is  indicated  as  D.  Since  the  signal  strength  of  the  reflections  is  sufficiently  small. 


27 


nonlinearities  are  not  reintroduced  on  the  return  path  through  the  medium.  The  RF  signal  present 
at  the  face  of  the  sensing  transducer  is  the  coherent  summation  of  each  harmonic  that  propagates 
back  to  the  sensing  transducer  after  signal  interaction  with  the  object.  The  RF  signal  is  illustrated 
in  D.  The  transducer  will  introduce  its  spatio-temporal  function  as  the  return  signal  is  sensed. 
This  reintroduction  of  transducer  effects  is  depicted  in  E.  In  E,  the  fundamental  and  second 
harmonic  signals  are  separated  using  simple  frequency  filtering  techniques.  The  fundamental 
signal  is  envelope  detected  at  the  fimdamental  frequency.  The  second  harmonic  signal  is  detected 
at  the  second  harmonic  frequency.  The  images  are  then  constructed  as  intensity  maps  of  the 
envelope  detected  signal. 

This  modular  approach  to  NUPROP  processes  allows  the  user  to  investigate  particular 
functional  dependences  relative  to  the  signal  at  various  places  within  the  propagation  path. 


Post  Processing:  Statistical  Analyses 

Statistical  analyses  are  perfijrmed  on  B-scan  unages  by  using  IDL  s  internal  function, 
MOMENT.  The  MOMENT  fimction  computes  the  mean,  variance,  and  skewness  of  a  sample 


population  contained  in  an  n-element  vector  X.  If  the  vector  contains  n  identical  elements. 


MOMENT  computes  the  mean  and  variance,  and  returns  the  IEEE  value  NaN  (not  a  number)  for 
the  skewness,  which  is  not  defined.  When  x  =  (xo,  Xj,  X2,  ....  xn.i),  the  various  moments  are 
defined  as  follows: 


_  1  AT-l 

MEAN  =  x  =  —Yx, 


(3-19) 


1  Jv-i  _ 

VARIANCE  =  =  -7— y!(^/  -  xf 

N-\1^ 


(3-20) 
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SKEWNESS  =  S^  =—y  - - 

Ni^[  o- 

An  image  is  converted  from  an  array  of  values  with  double  indices,  i  and  j,  to  its  lexicographic 
representation,  x*,  where  k  ranges  from  0  to  A  - 1  =  i  •  y  - 1 .  is  the  total  number  of  elements  in 

the  array  (or  vector). 

Along  with  the  statistics  of  the  fundamental,  second  harmonics,  and  third  harmonic  images, 
their  respective  intensity  histograms,  and  thresholded  images  are  presented.  The  histograms  give 
a  graphical  representation  of  the  statistical  values  associated  with  the  images.  Simple  thresholding 
of  the  images  show  the  potential  for  object  detection  and  segmentation. 


(3-21) 


simulated  images  for  the  fiindamental,  second  and  third  harmonics,  and  the  phantom  along  with  their  respective  thresholded 
images.  For  this  example,  a  simple  threshold  level  of  40  was  used.  The  upper  right  panel  displays  the  histogram  plots  for  each 
simulation  image.  The  lower  right  panel  displays  each  simulated  image’s  statistics:  mean,  variance,  and  skew. 
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Chapter  3:  Summary 

NUPROP  is  a  fairly  complex  simulation  and  analysis  tool  for  ultrasoxmd  modeling.  NUPROP 
represents  an  end-to-end  modeling  environment  fisr  studying  ultrasound  wave  propagation  and, 
ultimately,  B-Scan  image  formation.  This  chapter  discussed  the  generation  of  objects  called 
phantoms  using  Bamber’s  continuous  point  scatterer  concept.  The  tissue  characteristic  of  interest 
was  acoustic  impedance.  Input  signal  waveforms  used  in  NUPROP  are  either  sinusoidal  or 
frequency  modulated  sinusoidal,  either  apodized  (for  image  formation)  or  non-apodized  (single 
frequency  or  continuous  used  for  code  verification  and  harmonic  signal  analyses). 

Nonlinear  propagation  was  modeled  using  either  the  Burgers’  equation  or  Fubini  solution 
approaches.  Blackstock  developed  a  closed  form  solution  to  the  Burgers’  equation  that 
possessed  the  shock  distance  parameter,  I,  which  is  calculated  from  system  parameters  of 
velocity,  initial  wave  amplitude,  nonlinear  parameter  (fi),  and  fiindamental  frequency.  Likewise, 
the  Fubini  approach  possessed  this  same  shock  distance  limitation. 

Linear  propagation  was  modeled  using  Angular  Spectrum  and  Lommel.  Angular  spectrum 
has  the  advantages  of  being  applicable  throughout  the  propagation  path,  accommodate  non- 
radiaUy  symmetric  transducer  geometries,  and  accurately  computes  nonlinear  or  harmonic 
generation.  Angular  spectrum  requires  a  trade-off  between  computational  speed/volume  and 
accuracy/precision  in  computing  trans-axial  beam  profiles.  Lommel,  on  the  other  hand,  computes 
the  on-axis  longitudinal  beam  profile  and  the  trans-axial  beam  profiles  at  depth  with  accuracy  and 
precision.  Unfortunately,  Lommel  formulation  suffers  from  two  major  limitations:  radial 
symmetry  only  and  linear  propagation  only.  Lommel  does  not  predict  harmonic  signal  formation 
nor  does  it  accurately  compute  off-axis  beam  profiles  for  non-radiaUy  symmetric  geometries. 
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Lommel  does  not  accurately  compute  the  pressure  field  just  in  fi-ont  of  the  transducer.  A 
minimum  propagation  distance  is  required  before  Lommel  can  be  used. 

Image  generation  is  accomplished  using  the  Hilbert  transform,  which  converts  a  RF  signal  to 
an  “envelope”  detected  intensity  signal. 

Finally,  NUPROP  performs  some  rudimentary  post  processing  such  as  statistical  analyses, 
histogram  plotting,  and  simple  image  thresholding  based  on  the  statistics  of  the  images. 
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Chapter  4.  Code  Development  and  Description 

This  chapter  discusses  NUPROP’s  code  development,  the  use  of  Interactive  Data  Language 
Environment  for  developing  the  graphical  user  interfece;  describing  menu  operations,  radio  button 
operations,  or  parameter  input,  with  computational  considerations  where  appropriate. 

NUPROP  Hierarchy 

The  primary  tool  for  modeling  wave  propagation  in  biological  media  is  NUPROP.  NUPROP 
is  a  code  developed  using  Interactive  Data  Language  (IDL),  a  high-level  language  environment 
from  Research  Systems,  Inc.  IDL  is  a  complete  computing  environment  for  the  interactive 
analysis  and  visualization  of  data.  IDL  integrates  a  powerful,  array-oriented  language  with 
numerous  mathematical  analysis  and  graphical  display  techmques. 

NUPROP  is  widget-based  in  that  a  graphical  user  interface  (GUI)  provides  the  interactive 
platform  for  the  user  to  exercise  the  physics  associated  with  nonlinear  and  linear  wave 
propagation.  Widgets  (or  controls,  in  the  terminology  of  some  development  environments)  are 
simple  graphical  objects  such  as  pushbuttons  or  sliders  that  allow  user  interaction  via  a  pointing 
device  (usually  a  mouse)  and  a  keyboard.  IDL  graphical  user  interfaces  are  constructed  by 
combining  widgets  in  a  tree-like  hierarchy.  Each  widget  has  one  parent  widget  and  zero  or  more 
child  widgets.  There  is  one  exception:  the  topmost  widget  (called  a  top-level  base)  is  always  a 
base  widget  and  has  no  parent. 

Figure  4-1  shows  the  tree-like  hierarchy  of  NUPROP.  NUPROP  consists  of  four,  first-level 
bases:  Menu,  Button,  Parameter,  and  Display. 
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^  See  Figure  4-2  for  more  detail 
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Figure  4-1.  NUPROP  hierarchy  of  processing.  Each  first-level  base  is  designated  with  an  underline:  Menu,  Radio  Button, 
Parameter  Input,  and  Display.  Each  second-level  base  is  designated  by  the  origination  of  an  arrow  such  as  File  under  the  Menu 
Base  Operations.  Each  third-level  base  is  designated  by  smaller  print  text  such  as  Sinusoid  Chirp  of  Radio  Button  Base  as  well 
as  a  termination  and  origination  of  an  arrow.  Each  fourth-level  base  is  designated  by  small  print  italicized  such  as  Plot  Images 
under  Menu  Base/File/Save. 

The  Menu  Base  Operations  consists  of  three  second-level  bases:  File,  Execute,  and  Help.  File 
consists  of  three  third-level  bases:  Save,  Color,  and  Quit.  Save  consists  of  two  fourth-level  bases: 
plots  and  images.  Execute  is  described  later  in  this  section.  Help  invokes  a  Windows-based  help 
menu  describing  NUPROP  and  its  various  operations. 

The  Radio  Button  Operations  Base  consists  of  five  second-level  bases:  Input  Type,  Nonlinear 
Process  Type,  Linear  Process  Type,  Aperture  Type,  and  Shock  Limited  Processing.  Input  Type 
consists  of  two  third-level  bases:  sinusoid  and  chirp.  Nonlinear  Process  Type  consists  of  two 
third-level  bases:  Burgers’  and  Fubini.  Diffraction  Type  consists  of  two  third-level  bases:  Angular 
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Spectrum  and  Lommel.  Aperture  Type  consists  of  two  third-level  bases:  Circular  and 
Rectangular.  Shock  Limited  Processing  consists  of  two  third-level  bases:  Yes  and  No. 

The  Parameter  Input  Base  consists  of  thirty  second-level  bases: 

1 .  Time:  time  scale  factor  over  which  to  display  the  temporal  waveform  (in  seconds) 

2.  Distance:  maximum  ultrasound  wave  penetration  depth  (in  meters) 

3.  Velocity:  group  velocity  of  the  ultrasound  wave  (in  meters  per  second) 

4.  Transmit  Frequency:  oscillation  frequency  of  the  ultrasound  (in  Hertz) 

5.  Attenuation  Factor:  media  attenuation  (in  dB  per  cm  per  MHz) 

6.  Amplitude:  initial  ultrasoimd  wave  amplitude  (in  MPa) 

7.  B/A  Ratio:  nonlinear  parameter  (in  m^ s/kg) 

8.  #  of  Bessels:  number  of  terms  to  include  in  the  determination  of  nonlinearity 

9.  Window  Factor:  temporal  apodization  of  waveform 

10.  Size  (X-Y):  cross-track  resolution  (in  pixels) 

1 1 .  Size  (Z):  along-track  resolution  (in  pkels) 

12.  Z-Position:  observation  distance.  Actual  distance  is  calculated  by  multiplying  the  distance  by 
the  Z-position  and  dividing  by  the  along-track  resolution  (in  pixels) 

13.  Size  (t):  number  of  time  intervals  within  the  time  scale  fector 

14.  Beam  Profiles:  number  of  range  intervals  over  which  to  determine  beam  profiles 

15.  Aperture  Radius:  radius  of  the  transducer  (in  meters) 

16.  FOV:  field  of  view  associated  with  the  cross-track  resolution  (in  meters) 

17.  X-Size:  rectangular  aperture  size  in  the  x-direction  (in  meters) 

18.  Y-Size:  rectangular  aperture  size  in  the  y  -direction  (in  meters) 
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19.  Line  Trace:  used  for  sinusoid  input  only.  Determines  which  frequency  slices  to  display  in  the 
range  versus  frequency  plots.  Also  determines  which  range  profiles  to  use  when  calculating 
the  propagation  distance  influence  on  harmonic  amplitudes  (in  pixel  value) 

20.  Center  Frequency:  center  frequency  of  the  receiver  transducer  (usually  the  same  transducer  as 
the  transmit  transducer,  in  Hz) 

21.  Bandwidth:  bandwidth  over  which  the  transducer  responds  (in  Hz) 

22.  Transducer  Exponent:  the  transducer  response  is  Gaussian.  When  this  exponent  is  eqiial  to 
two,  then  the  transducer  is  the  typical  Gaussian  profile.  When  the  value  is  three  or  greater, 
the  Gaussian  response  changes  progressively  more  toward  a  RECT  fiinction. 

23.  Lower  Threshold:  Coupled  vdth  Upper  Threshold  establishes  the  range  over  which  pixel 
values  are  either  turned  black  or  white.  Used  in  statistical  analyses  to  determine  ease  of 
detection  based  on  differential  acoustic  impedance 

24.  Signal  to  Noise  Ratio:  Value  associated  with  the  power  ratio  of  signal  to  noise 

25.  Chirp  Start  Frequency:  start  frequency  of  the  chirp  input  signal  (in  Hz) 

26.  Chirp  Frequency  Modifier:  value  that  modifies  (increase  or  decreases)  the  starting  chirp 
frequency 

27.  Chirp  Gauss  Modifier  1:  the  chirp  input  signal  is  apodized.  In  combination  with  the  second 
Chirp  Gauss  Modifier,  this  Chirp  Gauss  Modifier  determines  the  overall  shape  of  the  Gaussian 
apodization  function 

28.  Chirp  Gauss  Modifier  2:  See  Chirp  Gauss  Modifier  1 

29.  Chirp  Time  Scale  Factor:  determines  the  temporal  vmdow  over  which  the  chirp  signal  is 
displayed 

30.  Chirp  Rate  Modifier:  modifies  the  chirp  rate. 
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The  Display  Base  consists  of  five  second-level  bases:  Temporal  Plot,  Frequency  Plot,  Range 
Images  (X-Z),  Profiles,  and  Angular  Image  (X-Y).  Range  Image  (X-Z)  consists  of  four  third- 
level  bases:  Fundamental,  Second  Flarmonic,  Third  Harmonic,  and  Phantom.  Profiles  consists  of 
five  possible  third-level  bases:  Fundamental,  Second  Harmonic,  Third  Harmonic,  Phantom,  and 
Beam. 

Figure  4-2  shows  the  continuation  of  the  Execute  option  under  Menu.  Execute  consists  of  six 
third-level  bases:  phantom  generation,  waveform  generation,  image  generation,  image  profiles, 
movies,  and  statistical  analyses.  Image  profiles  consists  of  three  fourth-level  bases:  x-z  image 
profile,  x-y  image  profile,  and  aperture  or  beam  profiles  (X-Y  Images).  Movies  option  consists  of 
two  fourth-level  bases:  temporal  movies  and  spatial  movies.  Temporal  movies  option  consists  of 
two  fifth-level  bases:  time  and  fi'equency  movies.  Spatial  movies  option  consists  of  two  fifth-level 
bases:  beam  profile  and  image  movies.  Image  movies  option  is  not  currently  active. 


Exectite  Operations 


Figure  4-2.  The  Execute  Base.  Execute  is  the  “heart”  of  NUPROP.  Through  the  Execute  Base  all  NUPROP  processes  are 
activated.  All  the  parameters  set  in  the  Parameters  Base  are  used  to  generate  the  waveforms  necessary  to  simulate  nonlinear 
ultrasound  propagation  and,  consequently,  image  formation. 
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Menu  Operations 

File 

Under  FUes  Operations,  the  user  can  implement  the  Save  routines  for  either  the  plots  or 
images  generated  from  the  Execute  operations.  Plots  consist  of  the  temporal  signals  (whether 
time  displayed  or  frequency  displayed).  The  angular  beam  profiles  are  also  saved,  if  generated 
during  the  Execute  operations.  Color  option  allows  the  user  to  choose  one  of  over  41  different 
color  schemes  built  into  IDL.  The  color  option  also  aUows  the  user  to  change  the  color  table  by 
introducing  a  user  defined  color  mapping.  The  Quit  option  is  self-explanatory. 

Execute 

The  user  can  generate  complex  phantoms  using  the  phantom  generate  option  under  execute. 
This  operation  initiates  a  module  called  NUPROP_Phantom_Generate.  This  module  allows  the 
user  to  modify  settings  for  three  elliptical  objects:  nominal  acoustic  impedance,  variation  in 
acoustic  impedance,  size,  shape,  and  orientation  of  the  objects.  The  user  can  change  the  variation 
in  acoustic  impedance  for  the  background  as  well.  The  simulation  phantom  and  its  pnstme 
representation  are  displayed.  The  user  can  save  the  simulation  phantom  and  its  pristine 
representation  by  selecting  the  save  option  under  file.  The  user  can  terminate  this  module  by 
selecting  quit  from  the  file  menu.  To  simulate  tissue,  digital  representations  are  generated  using 
Bamber’s  concept  of  continuous  array  of  point  sources  representing  the  phenomenon  of  mterest. 
NUPROP  simulates  acoustic  impedance  by  generating  a  complex  image  of  point  sources  as 

follows: 

1 .  Determine  the  number  and  size  of  scatter  regions 

2.  Determine  the  strength  of  acoustic  impedance  for  each  scatter  region 
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3.  Determine  each  region’s  variation  in  acoustic  impedance 

4.  Generate  the  image  based  on  these  determinations  using  an  IDL  code  called 
NUPROP_Tissue_Generate.  NUPROP_Tissue_Generate  is  a  module  incorporated  into 

NUPROP. 

Figure  4-3  shows  an  example  of  a  complex  phantom  generated  by  NUPROP.  NUPROP 
generates  three  elliptical  objects  based  on  the  parameters  set  by  the  user.  The  average  value 
of  the  tissue  simulated  image  is  zero.  This  prevents  the  introduction  of  an  artificially  induced 
acoustic  impedance  difference.  So,  the  variation  is  the  only  parameter  that  distinguishes  each 
scatter  region.  The  nominal  levels  are  used  to  generate  the  pristine  image  so  the  user  can 
visually  cue  the  anomalous  scatter  region(s).  The  major  and  minor  axis  modifiers  and  the 
object  extent  modifier  determine  the  overall  size  of  the  object.  The  object  rotation  angle  (in 
degrees),  object  shift  in  x,  and  object  shift  in  y  determine  the  positioning  and  orientation  of  the 
objects  in  the  256x128  array.  This  tissue  simulation  phantom  is  equivalent  to  Bamber’s 
continuous  array  of  point  scatterers. 


Three  elliptical  objects,  placed  and  oriented  in  a  256x128  array 


Figure  4-3.  Example  of  NUPROP’s  generation  of  tissue  phantoms.  Three  elliptical  objects  are  onented  and  inter-dispers^ 
within  a  256x128  array.  The  variation  determines  the  regional  distinction  in  the  tissue  simulation  (left  bottom  image).  e 
nominal  levels  determine  the  distinction  between  regions  in  the  pristine  simulation  (bottom  right  image).  The  pristine  image  is 
used  strictly  for  visually  cueing  of  anomalous  reglon(s).  This  example  shows  one  hypo-echoic  region  (cyst). 


The  waveform  execute  operation  begins  the  process  of  generating  nonlinear  waveforms. 
Input  from  the  Radio  Button  settings  and  parameter  values  determine  how  the  waveform  will  be 
generated.  Figure  4-4  illustrates  the  NUPROP  widget. 
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Figure  4-4.  Example  of  NUPROP  Widget-Based  Graphical  User  Interface.  This  example  illustrates  the  basic  widget 
configuration  including  the  radio  buttons,  parameter  input  section,  display  regions,  and  the  status  bar. 

The  radio  button  settings  determine  the  type  of  input  signal  (single  frequency  sinusoid  or 
chirp),  nonlinear  processing  (Burgers’  or  Fubini),  linear  or  diffraction  processing  (angular 
spectrum  or  Lommel),  aperture  type  (circular  or  rectangular),  and  whether  or  not  the  shock 
distance  limitation  will  be  invoked. 

The  parameter  settings  determine  the  fundamental  characteristics  of  the  output  signal.  The 
general  scenario  is  determined  by  the  time  and  distance  parameters.  The  material  type  is 
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determined  by  velocity,  B/A  Ratio,  and  attenuation  factor  parameters.  To  determine 
computational  size,  #  of  Bessels,  Size  X-Y,  Size  Z,  Size  t,  and  Beam  Profiles  parameters  are  used. 

Depending  on  the  radio  button  settings,  different  parameter  sets  within  the  parameter  widget 
are  used.  If  the  input  signal  is  a  single  fi-equency  sinusoid,  transmit  fi-equency,  wmdow  factor, 
center  fi-equency,  bandwidth,  transducer  exponent,  and  line  trace  parameters  are  used.  If  the 
input  signal  is  the  chirp,  then  chirp  start  fi-equency,  chirp  fi-equency  modifier,  chirp  Gauss 
modifiers  1  and  2,  chirp  time  scale,  chirp  rate  modifier,  and  transducer  parameters  are  used.  If  the 
transducer  geometry  is  circular,  then  the  aperture  radius  and  FOV  are  used.  If  the  transducer 
geometry  is  rectangular,  then  the  X-Size  and  Y-Size  are  used  along  with  the  FOV  parameter. 

The  user  can  generate  simulated  ultrasound  images  by  selecting  image  generate  under  the 
execute  option.  NUPROP  will  use  the  parameter  settings  that  generated  the  waveforms.  The 
simulated  ultrasound  images  will  be  displayed  within  the  RANGE  IMAGE  (X-Z)  region  of  the 
widget  (see  Figure  4-4)  for  the  fundamental,  second  and  third  harmonics,  and  the  phantom  used. 
The  transducer  geometry  is  displayed  within  the  ANGULAR  IMAGE  (X-Y)  region. 

Statistical  analyses  can  be  performed  on  the  images  by  selecting  Statistical  Analyses  m  the 
execute  menu.  A  display  of  the  statistics  (mean,  variance,  skew),  a  histogram  of  the  mtensity 
values,  and  a  rudimentary  threshold  segmentation  of  the  fundamental,  second  and  third  harmomcs, 
and  the  phantom  images  are  presented.  Figure  4-5  (repeat  of  Figure  3-4)  illustrates  an  example 

statistical  analyses. 
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Figure  4-5.  Example  of  statistical  ansilyses  from  NUPROP  ultrasound  simulation  widget.  The  left  panel  displays  the 
simulated  images  for  the  fundamental,  second  and  third  harmonics,  and  the  phantom  along  with  their  respective  thresholded 
images.  For  this  example,  a  simple  threshold  level  of  40  was  used.  The  upper  right  panel  displays  the  histogram  plots  for  each 
simulation  image.  The  lower  right  panel  displays  each  simulated  image’s  statistics:  mean,  variance,  and  skew. 


Computational  Considerations  Associated  with  Radio  Button  Selections 
Input  Type 

If  the  windowing  factor  is  sufiBciently  small,  on  the  order  of  0.00005,  then  a  single  sinusoid 
input  is  essential  continuous.  However,  if  the  windowing  factor  is  on  the  order  of  0.5,  then  the 
single  sinusoid  input  is  pulsed  (finite  temporal  duration).  A  pulsed,  single  sinusoid  input  is  no 
longer  monotonic  but  is  comprised  of  a  “band”  of  Ifequencies  centered  on  the  fimdamental 
fi-equency.  This  characteristic  significantly  impacts  the  computational  results:  onset  distances  for 
second  and  third  harmonics,  and  image  formation.  For  the  chirp,  if  the  windowing  factor  is  on  the 
order  of  5xl0■’^  then  the  input  signal  is  essentially  continuous.  If  the  windowing  factor  is  on  the 
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order  of  5x10 then  the  chirp  becomes  an  apodized  or  pulsed  chirp  which  translates  to  a  “low 
passed”  chirp  signal  in  the  frequency  domain. 

Nonlinear  Processim:  Burners '  or  Fubini 

Equation  3-4  gives  the  closed  form  solution  to  Burgers’  equation  as  derived  by  Blackstock. 
This  solution  incorporates  the  limitation  of  the  shock  distance  in  the  calculation  of  the  distorted 
waveform.  NUPROP  will  not  compute  the  waveform  unless  the  shock  distance  limitation  is 
turned  off.  Once  the  shock  distance  limitation  is  turned  off,  then  NUPROP  computes  the 
waveform.  The  user  should  realize  that  the  waveform  computed  is  no  longer  valid.  The  Fubim 
approach  also  possesses  this  shock  distance  limitation. 

Linear  Processing:  Lommel 

Computationally,  Lommel  is  significantly  fester  than  angular  spectrum.  However,  the  Lommel 
formulation  approach  is  strictly  a  linear  processing  technique  and  cannot  be  used  for  forward 
propagation  of  the  ultrasound  wave,  since  the  medium  introduces  nonlinearities,  provided  the 
scenario  parameters  are  conducive  to  harmonic  wave  formation.  To  correctly  predict  the 
waveform  in  the  outbound  direction,  a  step  by  step  (small  increments)  approach  is  required. 
Angular  spectrum  provides  this  incremental  determination  of  the  waveform.  The  field  at  a 
propagated  distance,  ,  depends  on  the  waveform  frequency  content  at  the  previous  step,  z„_, . 
Once  the  field  has  interacted  with  the  target  region,  the  return  signal  is  sufficiently  small  that 
nonlinearities  are  not  introduced  on  the  return  or  inbound  direction.  Lommel  can  be  used  for  the 
inbound  or  return  propagation.  Additionally,  the  Lommel  formulation  used  in  NUPROP  requires 
circular  (radial)  symmetry.  Under  very  restricted  conditions,  however,  the  Lommel  formulation 
can  produce  reasonable  approximations  for  on-axis  signal  amplitudes  for  rectangular  arrays. 
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Shock  Limited  Distance  Calculation 


As  mentioned  previously,  the  shock  distance  limits  the  range  over  which  valid  computations 
can  be  made.  NUPROP  allows  the  user  to  circumvent  this  limitation  with  the  understanding  that 
observation  distances  beyond  the  shock  distance  will  produce  invalid  results. 

Computational  considerations  associated  with  parameters 
Format  Size 

The  format  size  greatly  impacts  the  con^utational  speed  of  NUPROP.  NUPROP  generates 
volumetric  data:  x-y-z  and  t.  The  waveforms  generated  using  NUPROP’s  nonlinear  module  are 
dependent  on  z  and  t.  The  images  generated  using  NUPROP’s  image  generation  module(s)  are 
dependent  on  x,  y,  z,  and  t.  As  an  example,  if  the  user  chooses  an  angular  array  size  of 256  x  256 
and  a  number  of  beam  proffles  of  64,  then  a  256  x  256  x  64  “data  cube”  will  be  generated.  This 
data  cube  is  for  storing  the  necessary  information  for  several  other  operations  such  as  beam 
movies  and  angular  spectrum  displays.  A  data  cube  of  this  size,  given  8-bit  precision  results  in  a 
3.36  Mbyte  data  cube.  If  the  angular  array  size  is  doubled,  then  the  data  cube  size  becomes  13.42 
Mbytes.  Care  must  be  exercised  when  generating  large  computational  arrays.  Memory  allocation 
limitation  as  well  as  computational  time  must  be  considered. 

Line  Trace 

The  line  trace  feature  allows  the  user  to  pick  a  set  of  frequency  traces  from  the  nonlinear 
waveform  generated  data  set.  Line  trace  is  appUcable  to  single  frequency  input  only.  Usually,  the 
line  trace  is  an  integer  divisor  of  the  fundamental  frequency.  For  instance,  using  a  sinusoidal  input 
at  2  MHz,  1500  m/s  wave  velocity,  time  of  20  ps,  line  trace  to  use  is  20.  The  list  below  provides 
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the  appropriate  values  for  the  line  trace  given  the  input  frequency,  with  a  wave  velocity  of  1500 
m/s,  time  of  20  ps. 


Frequency  (MHz) 

Line  Trace  (Pixel  Number) 

0.8 

8 

1.0 

10 

2.0 

20 

3.0 

30 

4.0 

40 

The  user  has  the  option  of  picking  any  line  trace  with  the  foUovring  limitation:  the  highest 
harmonic  line  trace  must  be  less  than  the  array  size  aUotted  for  frequency  computations.  For 
example,  if  the  frequency  computational  array  is  set  to  512,  and  the  number  of  frequency  terms  to 
compute  is  three,  then  the  third  harmonic  cannot  exceed  512,  which,  in  turn,  limits  the  line  trace 

value  to  170. 

Number  of  Bessel  Terms. 

The  number  of  Bessel  terms  impacts  the  computation  of  the  nonlinear  waveform.  NUPROP 
uses  the  number  of  Bessel  terms  to  determine  the  number  of  frequency  terms  to  compute.  This 
process  is  relatively  fast,  but  is  platform  dependent.  For  a  350MHz  Pentium  II,  the  runtime  for 
computing  30  frequency  components  (each  component  waveform  vector  size:  512)  is 
approximately  35  seconds.  AdditionaUy,  the  number  of  Bessel  terms  computed  must  take  into 
account  the  temporal  frequency  array  size.  If  the  array  size  is  too  smaU,  then  higher  order  Bessel 
terms  will  be  aliased.  Typical  computational  size  settings  are  as  follows: 

1.  #  of  Bessels:  3  -  since,  for  most  harmonic  signal  generation  conditions,  90%  or  more  of  the 
wave  energy  is  captured  in  the  first  three  terms  -  fimdamental,  second  and  third  harmonics 

2.  Size  X-Y:  256  -  this  array  size  yields  reasonably  accurate  cross-track  beam  profiles 
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3.  Size  Z:  256  -  this  vector  size,  combined  with  Size  X-Y,  yields  a  reasonably  sized  volumetric 
array  for  accuracy,  computational  memory  requirements,  and  computational  time 

4.  Size  t:  512  -  this  vector  size  yields  reasonably  accurate  temporal  and  frequency  plots, 
especially  for  #  of  Bessels  parameter  size  in  excess  of  10 

5.  Beam  Profiles:  8  -  this  vector  size  is  suitable  for  “gross”  approximations  of  down-range,  on- 
axis  profiles.  For  better  accuracy,  a  Beam  Profiles  value  of  128  is  appropriate.  The 
computational  time  will  greatly  increase  (on  the  order  of  five  to  six  minutes  to  compute  the  beam 

profile). 

Chapter  4:  Summary 

In  summary,  NUPROP  is  a  modularly  design  IDL  code  that  models,  end-to-end,  nonlinear 
ultrasound  wave  propagation  and  ultrasound  image  formation.  NUPROP  is  based  on  Svidgets, 
or  controls,  in  the  terminology  of  some  development  environments.  Widgets  are  simple  graphical 
objects  such  as  pushbuttons  or  sliders  that  allow  user  interaction  via  a  pointing  device  (usually  a 
mouse)  and  a  keyboard.  IDL  graphical  user  interfaces  are  constructed  by  combining  widgets  in  a 
tree-like  hierarchy.  Each  widget  has  one  parent  widget  and  zero  or  more  child  widgets.  There  is 
one  exception:  the  topmost  widget  (called  a  top-level  base)  is  always  a  base  widget  and  has  no 
parent.  So,  NUPROP  is  a  graphical  user  interface  (GUI)  providing  the  interactive  platform  for 
the  user  to  exercise  the  physics  associated  with  nonlinear  and  linear  wave  propagation. 

NUPROP  is  built  with  modules  (operations  that  are,  essentially,  stand-alone  in  that  they  do 
not  need  the  GUI  to  work)  and  fiinctions  (modular  code  that  requires  parameter  values  passed 
through  the  GUI  to  operate).  The  complexity  of  each  module  or  function  determines  the 
computational  requirements,  hence  the  computational  time,  accuracy,  and  precision. 
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Chapters.  NUPROP Verification 

This  chapter  describes  the  verification  of  NUPROP  as  compared  to  existing  data,  other  code 
simulations,  and  conditions  required  to  generate  nonlinearly/linearly  propagated  ultrasound  waves. 
Sections  are  divided  into  Longitudinal  Signal  Formation  (nonlinear  propagation),  Trans-axial 
Signal  Formation  (linear  or  diffraction),  and  Image  Formation. 


Longitudinal  Signal  Formation. 

Figure  5-1  shows  the  comparison  between  NUPROP’s  Fubini  and  Burgers’  solutions  for  wave 
propagation  in  water  at  23°c,  at  5  centimeters  and  9.52  centimeters  fi-om  the  transducer  (2.58 
MHz,  circular)  and  data  fi-om  Ryan,  1963. 


Figure  5-1.  Comparison  of  NUPROP  output  and  data  from  Ryan,  R.  P.,  Ph.D.  thesis.  Brown  Uni  versit)^,  1963.  tte  un Jrl^  g 
image  is  data  captured  from  photographing  the  oscilloscope  output  at  two  penetration  depths;  5  and  10  cm. 
traSs  are  NUPROP  predictions  using  Fubini  (black  dotted)  and  Burgers’  (dotted  red)  under  the  same  (5  cm  penetration  depth) 
or  similar  conditions  (9.52  cm  vice  10  cm).  Note  that  under  conditions  prescribed  for  (b),  WROP  predicts  the  nonline^ 
waveform  precisely.  For  condition  (c),  NUPROP  predicts  the  waveform  at  the  shock  distance  (9.52  cm),  comparing  favorably 

with  Ryan’s  data  at  10  cm. 
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Ryan  defined  eras  the  ratio  of  the  observation  distance  to  the  shock  distance.  In  the  case  of 
condition  (c)  x  =  10  centimeters,  a  is  1.05,  slightly  passed  the  computed  shock  distance,  9.52 
centimeters.  NUPROP  checks  to  see  if  the  observation  distance  exceeds  the  shock  distance.  If  it 
does,  then  NUPROP  will  not  compute  the  waveform  since  the  model  is  invalid  beyond  the  shock 
distance.  However,  the  user  can  bypass  this  limitation  by  selecting  “no”  for  shock  limit  (radio 
button  option  in  NUPROP),  realizing  that  the  results  are  not  valid.  For  condition  (b),  NUPROP 
results  (Fubini  and  Burgers’)  are  nearly  an  overlay  of  the  experimental  data.  NUPROP  used  15 
fi-equency  components  to  generate  the  waveform.  According  to  Christopher  (personal 
communications,  1998),  no  less  than  20  to  25  terms  should  be  used  to  reduce  any  “wavy”  artifact 
in  the  reconstruction.  For  condition  (c),  NUPROP  produced  results  within  1%  error  using  Fubini 

or  Burgers’. 

Figure  5-2  shows  the  growth  and  decay  of  harmonics  for  2.5MHz  ultrasound  m  water.  The 
comparison  is  between  data  coUected  by  Ryan,  et.  al.  (left  graph)  and  results  from  simulations 
with  NUPROP  (right  graph).  Ryan’s  theoretical  predictions  are  the  solid  curves.  The  other  data 
in  Ryan’s  graph  are  (x)  digital  analysis;  (X)  analyzer  measurement.  Digital  analysis  was 
accomplished  via  photographing  waveforms  on  the  oscilloscope  and  digitizmg  the  photographs. 
Analyzer  measurements  were  accomplished  by  passing  detected  signal  through  electric  circuit 
tuned  to  the  appropriate  harmonic.  The  graph  on  the  right  is  from  NUPROP  under  the  same 
conditions  from  Ryan,  et.  al. :  2.5  MHz  ultrasound  wave  in  water.  NUPROP  settings  were  wave 
velocity:  1550  m/s;  number  of  Bessel  terms:  4;  B/A  ratio:  5.3.  NUPROP  predictions  are  less  than 
1  percent  deviation  from  Ryan’s  predicted  growth  and  decay  for  the  fundamental,  second,  and 

third  harmonics. 
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Figure  5-2.  Comparison  of  fundamental  decay  and  harmonic  growth.  Data  on  the  left  is  from  Ryan.  R.  P.,  A.  Lutsch,  and  R.  T. 
Beyer,  J.  Acoust.  Soc.  Am.,  34,  31,  (1962). 

Figure  5-3  shows  NUPROP’s  predicted  on-axis  longitudinal  signal  amplitude  (fundamental 
and  second  harmonic)  as  a  function  of  propagation  distance  in  comparison  to  data  from  Averkiuo, 
Roundhill,  and  Power,  1997.  Averkiuo,  et.  al.  assumed  that  a  circular  transducer  with  equivalent 
area  to  a  rectangular  transducer  would  produce  equivalent  on-axis,  longitudinal  signal  ampUtude 
profiles.  The  solid  black  line  represents  Averkiuo’s  2  MHz  fundamental  on-axis  longitudinal 
profile  for  a  P3-2  geometry.  The  thick  dashed  line  is  NUPROP’s  prediction  using  an  actual  P3-2 
transducer  configuration.  Likewise,  the  thin  dashed  line  is  Averkiuo’s  second  harmonic  on-axis 
longitudinal  profile.  The  thin  dotted  line  is  NUPROP’s  second  harmonic  axial  amplitude  profile 
using  a  P3-2  transducer  configuration. 
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Figure  5-3  Comparison  of  data  from  Averkiuo  et.  al.  and  NUPROP  for  a  P3-2  rectangular  geometry  transducer  on-axis  signal 
strength  as  a  function  of  propagation  distance.  The  solid  line  represents  Averkiuo’s  on-axis  signal  amplitude  at  the 
fundamental  frequency  of  2  MHz,  10  cm  focal  length,  circular  transducer  with  equivalent  area  to  the  P3-2  rectangular 

transducer. 

Averkiuo  assumed  a  circular  aperture  of  equivalent  area  would  produce  an  equivalent  on-axis 
signal  attenuation  as  an  actual  P3-2  geometry.  The  difference  in  second  harmonic  traces  is  due  to 
the  differences  in  using  circular  equivalence  and  an  actual  P3-2  configuration.  For  the 
fimdamental,  the  difference  is  minimal,  whereas  the  difference  in  second  harmonic  is  more 
significant,  ranging  upwards  of  20  percent.  The  “wavy”  nature  of  the  fundamental  for  NUPROP 
is  due  largely  to  wrap-around  effects  (or  spatial  errors).  It  should  be  noted  that  the  spatial 
wavefi'ont  is  not  attenuated  at  the  edges,  as  in  both  Averkiuo  s  and  Christopher  s 
implementations.  The  spatial  apodization  will  eliminate  the  “wavy”  nature  of  the  fundamental  as 
well  as  correct  for  the  differences  shown  for  the  second  harmonic  in  Figure  5-3.  This  data 
supports  Averkiuo’s  supposition  of  equivalence,  but  does  reveal  potentially  large  errors, 
especially  for  second  harmonic  amplitudes. 
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Figure  5-4  further  reveals  the  danger  of  Averkiuo’s  equivalence  approximation.  As  aperture 
size  increases,  errors  slightly  decrease:  8  pixels  diameter  produced  a  maximum  error  of 
23.0469%,  16  pixels  produced  a  maximum  error  of  22.1680%,  and  32  pixels  produced  a 
maximum  error  of 21.6553%. 


Figure  5-4.  Demonstration  of  difference  between  circular  and  rectangular  aperture  diffraction  patterns  and  its  effect  on  signal 
strength.  As  the  size  of  the  aperture  increases,  the  maximum  signal  difference,  in  percent,  slightly  decreases.  8  pixels  produced 
a  maximum  difference  of  23.0469%;  16  pixels  produced  a  maximum  difference  of  22.1680%;  32  pixels  produced  21.6553% 
difference. 

Figure  5-5a,b  shows  a  comparison  between  NUPROP’s  angular  spectrum  and  Lommel 
formulation  longitudinal  on-axis  signal  attenuation  as  a  fimetion  of  propagation  distance.  Angular 
spectrum  technique  has  a  computational  limitation  that  requires  a  minimum  angular  resolution  of 
256  by  256  pixels.  Figure  5-5a  illustrates  the  effect  on  increasing  the  angular  resolution  from  32 
pixels  to  256  pixels,  keeping  the  range  resolution  fixed  at  8  pixels.  The  solid  black  curve  is  the 
Lommel  formulation  solution  for  on-axis  signal  strength  as  a  fimetion  of  propagation  distance.  A 
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circular  piston  transducer  with  1.0  cm  diameter,  center  frequency  0.956  MHz,  range  of  23.9  cm, 
and  the  wave  group  velocity  of  1200  m/s  was  used.  At  coarse  range  resolution  (8  datapoints  over 
the  entire  range),  the  angular  spectrum  techmque  appears  to  “follow  the  average  value  of 
Lommel  formulation  quite  well  (256  x  256  cross-track  resolution).  It  should  be  noted  that  the 
Lommel  curve  used  1024  datapoints  in  the  down-range  dimension  (z-direction).  However,  as  the 
cross-track  resolution  decreases,  the  angular  spectrum  representation  begins  to  deviate 
substantially.  Figure  5-5b  shows  the  effect  of  increasing  down-range  resolution  while  maintaining 
the  cross-range  resolution  fixed.  As  is  readily  apparent,  the  down-range  resolution  has  little  effect 
on  the  overall  accuracy  of  the  angular  spectrum.  Cross-track  resolution  significantly  impacts  the 
fidelity  of  the  representation  relative  to  the  Lommel  formulation  values. 

(a)  (t>) 


Figure  5-5  Comparison  of  angular  spectrum  (Red  Curves)  and  Lommel  (Black  Curve)  formulation  accuracy,  (a)  Effects  of 
increasing  cross-track  resolution  on  angular  spectrum  accuracy,  (b)  Effects  of  down-range  resolution  on  angular  spectrum 
accuracy. 
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Trans-Axial  Signal  Formation 

Trans-axial  signals  are  formed,  primarily,  by  diffraction.  Figure  5-6  compares  NUPROP’s 
generation  of  trans-axial  beam  profiles  and  results  fi-om  Christopher’s  NLP  code.  The  black 
curves  correspond  to  Christopher’s  NLP  code  and  the  red  curves  correspond  to  NUPROP’s 
Fubini- Angular  Spectrum  Techmcjue.  For  trans-axial  distances  less  than  1.5  cm,  NUPROP  and 
NLP  are  within  a  few  percent.  Beyond  1.5  cm,  NUPROP  and  NLP  results  differ  as  much  as  50% 
in  amplitude  value.  The  side  lobe  termination  points  (valley  points)  are  the  same  for  both 
NUPROP  and  NLP.  The  differences  in  anplitude  values  for  trans-axial  distances  greater  than  1.5 
cm  is  due  to  the  finite  extent  in  the  NUPROP  calculations.  NUPROP  assumes  a  finite  boundary, 
or  distance,  over  which  to  calculate  the  trans-axial  beam  profile.  This  boundary  is  manifested  by 
the  field  of  view  parameter  causing  a  “wall  effect”  or  wrap-around  to  occur,  thus  altering  the 
overall  amplitude  of  beam  profile  outside  of  1 .5  cm.  This  difference  is  more  prominent  in  the 
second  harmonic  beam  profile  than  the  fundamental.  These  differences  can  be  corrected  by  using 
a  spatial  apodization  to  smooth  out  the  edges,  thus  eliminating  the  wrap-around  in  the 
calculations. 


52 


NLP  -  Black  Lines  NUPROP  -  Red  Lines 


Trans-axial  Distance  (m) 


Figure  5-6.  Comparison  of  NUPROP’s  Fubini-Angular  Spectrum  technique  with  Christopher’s  NLP  Code  for  trans-axial  signal 
formation  (Beam  Profile).  Conditions:  2  MHz  fundamental  with  harmonic  at  4  MHz,  P3-2  transducer  in  water,  at  4  cm 
penetration.  The  black  curves  were  generated  by  Christopher’s  NLP  code.  The  red  curves  were  generated  using  NUPROP’s 
Fubini-Angular  Spectrum  Technique. 

Another  feature  readily  apparent  from  Figure  5-6  is  the  trans-axial  distance  associated  with 
the  main  lobes  of  the  fundamental  and  harmonic  beam  profiles.  The  second  harmonic  main  lobe’s 
trans-axial  extent  is  approximately  half  of  the  fundamental,  suggesting  that  the  second  harmonic 
will  yield  improved  lateral  resolution.  It  should  be  noted  that  the  second  harmonic  is  normalized 
relative  to  its  maximum  value,  readily  comparing  the  main  lobe  beam  widths.  Typically,  the 
second  harmonic  signal  strength  is  several  dB’s  lower  than  the  fundamental. 


Image  Formation 

Figure  5-7  shows  a  comparison  of  NUPROP’s  image  generation  for  a  point  source  under  the 
conditions  specified  in  Bamber  and  Dickinson  [1980].  The  images  on  the  left  are  for  envelope 
detected  signal  and  RF  signal  (top  and  bottom,  respectively).  The  images  on  the  right  are  from 
NUPROP’s  simulation  of  point  source  images,  envelope  detected  and  RF  (top  and  bottom. 
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respectively).  The  white  boundaries  designate  the  3dB,  or  half  power,  contour  for  the  envelope 
and  RF  signals,  assuming  symmetric  pulse  shape.  Bamber  and  Dickinson’s  data  exhibit  a 
“skewed”  distribution  because  of  the  non-symmetric  rise  and  fall  time  of  their  temporal 
apodization  function.  The  inscribed  contour  is  associated  with  the  more  gradual  fall  off. 
NUPROP  2issumes  a  symmetric  rise  and  fall  time  for  its  temporal  apodization  fiinction.  As  is 
readily  apparent  by  this  comparison,  NUPROP  produced  RF  and  envelope  detected  signals 
comparable  with  Bamber  and  Dickinson. 

Date  foan  Bamber  and  Data  from  NUPRDP  Using 

Diddnson  (1 980)  Bamber’s  Conditioiis 

Envdope  Detected  Signal 


RF  Signal  for  Point  Source 


Figure  5-7.  Comparison  of  NUPROP  simulation  of  a  point  source  signal  using  conditions  described  in  Bamber  and  Dickinson 
[1980]  and  data  from  Bamber  and  Dickinson  [1980].  The  inscribed  boundaries  (white  circle  [top]  and  ellipse  [bottom])  are 
associated  with  the  3dB,  or  half  power,  contours. 

Chapter  5:  Summary 

NUPROP’ s  nonlinear  wave  generation  and  propagation  results  compared  favorably  with  data 
found  in  Uterature,  having  less  than  one  percent  difference  comparing  NUPROP’s  Fubini 
implementation  with  existing  data  (Ryan,  [1963]  and  less  than  a  few  percent  between  NUPROP’s 
Burgers’  solution  and  existing  data.  NUPROP  also  accurately  predicted  on-axis  wave  amplitude 
for  propagation  distances  comparable  to  Averkiuo’s  simulation  conditions.  Additionally, 
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NUPROP  confirmed  Averkiuo’s  supposition  that,  for  equivalent  areas,  circular  transducers 
produce  equivalent  on-axis  longitudinal  amplitudes  as  rectangular  arrays  (P3-2  geometry). 
However,  NUPROP  also  showed  that  for  wave  or  field  amplitudes  off-axis,  Averkiuo’s 
assumption  is  invalid. 

In  comparing  trans-axial  beam  profiles  derived  fi'om  diffraction  coupled  with  nonlinear 
propagation,  NUPROP  compared  very  favorably  with  literature  (Christopher,  [1997]),  having  less 
than  a  few  percent  difference  in  trans-axial  beam  profile  amplitudes  fiir  2  MHz  and  its  second 
harmonic  at  4  MHz,  within  1.5  cm  radial  distance  fi-om  the  axis. 

Comparing  NUPROP’s  RF  and  envelope  detection  simulated  images  with  Bamber’s  teehnique 
for  a  point  source.  Bamber’s  simulation  used  a  “skewed”  or  non-symmetric  temporal  pulse  shape, 
whereas  NUPROP  assumes  a  temporally  symmetric  beam  apodization.  Comparing  the  3  dB  point 
for  the  most  gradual  M-off,  NUPROP  is  within  one  percent  in  computing  the  half  power  point, 
or  3  dB  down  fi-om  maximum  intensity. 

Given  the  favorable  comparison  of  NUPROP’s  simulations  of  nonlinear  wave  generation  and 
propagation  using  Fubini  or  Burgers’  partial  differential  equation  solutions,  linear  wave 
propagation  (diffraction)  using  angular  spectrum  or  Lommel  formulation,  and  image  formation 
using  Bamber’s  technique,  verification  of  NUPROP’s  viability  as  a  simulation  code  is  complete. 
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Chapter  6.  Exercising  NUPROP 

This  chapter  discusses  variations  of  NUPROP  parameters  and  the  effect  of  variations  on  signal 
and  image  generation.  Parameters  of  interest  include  computational,  scenario  or  environmental, 
object,  and  system  parameters  as  depleted  in  Figure  6-1 . 


Ccmputational: 

-  Grid  Size(x,y,z:t) 

-  #af  Bessel  Terms 

-  #  of  Profiles 

-  Statistical  Analyses 


Eiivironmental: 

-  Homogeneous  Path 

-  B/A  Ratio 

-  Attenuation 


Syston:  Object; 

-  Chirp  or  Sinusoid  .  b/A  Ratio 

-  Burger’s  vs.  Fubini  .  Acoustic  In^edance 

-  Transducer  Characteristics 

-  Frequaicy 

-  An5)litude 

-  Pulse  Shaping 

-  Observation  Distance 

-  Frequency  Trace  Lines 


Figure  6-1.  Graphical  representation  of  NUPROP  variables.  Computational  variables  govern  the  size  of  the  computational 
space  allocated  for  signal  and  image  generation.  Environmental  variables  govern  path  constituenQ-  and  attenuation.  Object 
variables  govern  signal  interaction.  System  variables  govern  input  type,  nonlinear  and  linear  propagation  modeling  types, 
transducer  characteristics,  path  observation  depth,  and  trace  lines  for  frequency  representational  graphs.  Image  of  heart  from 
Gray’s  Anatomy.  New  York:  Bounty  Books,  p.  467,  1977. 

For  signal  generation  results,  namely  onset  distance  computations,  no  longitudinal  or  temporal 
apodization  effects  are  imposed.  Thus,  the  input  signal  is  a  single  frequency  sinusoid.  For  image 
formation  results,  a  pulsed  or  narrowband  sinusoid  input  signal  is  used.  Care  is  taken  to  ensure 
that  the  narrowband  results  do  not  possess  frequency  overlapped  signals. 
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By  this,  the  second  and  subsequent  harmonic  signals  are  band  limited  in  such  a  way  as  to 
reduce  the  potential  of  transferring  energy  between  frequencies.  Figure  6-2  illustrates  this 
concept.  Assuming  that  the  temporal  apodization  is  Gaussian  with  exponent  2,  then  the  band 
about  each  frequency  must  be  constructed  such  that  the  overlap  between  frequency  bands  is 
minimized.  This  will  ensure  that  the  assumption  of  coherent  addition  applies  to  the  sum  signal  for 
the  overall  signal  comprised  of  fundamental,  second,  and  subsequent  harmonics. 


Third  Harmonic 


Figure  6-2.  Dlustration  of  minimizing  frequency  band  overlap  to  ensure  minimum  energy  cross-over  between  frequencies.  By 
minimizing  the  overlaps,  the  coherent  summation  of  frequency  components  within  the  nonlinearly  generated  signal  is  valid. 

Single  Frequency  Sinusoid:  Onset  Distance  Computation 

This  section  discusses  the  effects  of  varying  B/A  ratio  and  path  attenuation  on  signal  and 
image  generation.  For  examples  in  this  section,  the  following  parameters  were  held  eonstant  and 
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constituted  a  point  of  reference  for  signals  and  images  generated  varying  B/A  ratio  and  path 


attenuation: 


Parameter 

Value 

Time: 

40  ps 

Distance: 

10  cm 

Velocity: 

1500  m/s 

Transmit  Frequency: 

1.0  MHz 

Amplitude: 

0.5  MPa 

#  of  Bessel  Terms: 

4 

Window  Factor: 

0.5 

XY  Size: 

256 

Z  Size: 

256 

Observation  Distance: 

32 

T  Size: 

512 

#  of  Beam  Profiles: 

8 

Aperture  Radius: 

6.5  mm 

FOV: 

10  cm 

X-Size: 

22.5  cm 

Y-Size: 

15  cm 

Line  Trace: 

50 

Transducer 

Center  Frequency: 

2.5  MHz 

Bandwidth: 

1.5  MHz 

Exponent: 

2 

SNR 

5x10’ 

Chirp 

Start  Frequency 

0.5  MHz 

Gauss  Modifier 

3.5x10'' 

Time  Scale  Factor 

1.0 

Rate  Modifier 

1.5x10" 

Continuous  Wave 

Effects  of  B/A  Ratio  and  Path  Attenuation  Variation  on  Onset  Distances  for  Second  and 
Third  Harmonic  Signal  Generation 

One  measure  of  harmonie  frequeney  generation  is  the  onset  distance.  Onset  distance  is 
defined  as  the  distance  at  which  a  given  harmonic  signal  first  reaches  five  percent  of  its  maximum 
value.  Figures  6-3  through  6-5  are  log-linear  plots.  If  changes  in  parameter  values  produce  linear 
dependence  of  onset  distance  in  this  representation  space,  then  the  process  is  said  to  be  first  order, 
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with  the  slope  of  the  trace  being  a  measure  of  the  rate  constant.  First  order  processes  can  be 
represented  as  a  negative  exponential 


Taking  the  natural  logarithm  of  both  sides  yields 


(6-1) 


(6-2) 


where  Pc  is  the  process  rate  constant.  So,  graphing  the  logarithm  oif(P)  as  a  function  of  ^yields 
plots  represented  in  figures  6-3  through  6-5,  where  the  slope  is  inversely  proportional  to  the  rate 
constant  associated  with  the  process. 

Figure  6-3  shows  the  effect  of  varying  B/A  ratio  on  the  onset  distance  for  a  pulsed  sinusoid, 
IMHz  center  fi-equency,  10  ps  pulse  duration.  Path  attenuation  was  set  to  1x10'®  dB/cm/MHz. 


B/A  ratio  was  varied  fi'om  5.0  to  12.0. 


_o_2nd 

-x-3rd 


B/A  Ratio 


Figure  6-3.  Effect  of  varying  B/A  ratio  on  harmonic  signal  onset  distance.  B/A  ratios  range  from  5.0  to  12.0,  typical  range  of 
biological  medium  (distilled  water  to  fat).  Onset  distances  range  from  4.7  mm  to  3.9  mm  for  second  harmonic  signal.  Onset 
distances  range  from  2.1  cm  to  1.7  cm  for  third  harmonic  signal. 
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These  results  suggest  that  onset  distance  for  second  and  third  harmonics  of  1  MHz 
fundamental  insonification  are  not  strongly  dependent  on  B/A  ratio  for  typical  biological  medium 

t 

nonlinear  parameter  values. 

Figure  6-4  shows  the  effect  of  varying  path  attenuation  from  IxlO'^  to  2.0  dB/cm/MHz  and 
B/A  ratio  from  5.0  to  12.0. 


B/A  Ratio 


Figure  6-4.  Effect  of  varying  B/A  ratio  from  5.0  to  12.0,  for  varying  path  attenuation,  on  onset  distance  for  pulsed  sinusoid 
with  center  frequency  of  1.0  MHz,  pulse  duration  of  10  ps. 

These  results  suggest  that  as  the  path  attenuation  increases,  the  onset  distances  decrease,  but 
only  after  attenuation  reaches  or  exceed  unity.  Though  the  curves  for  each  path  attenuation 
decreases  as  B/A  ratio  increases,  the  slope  is  relatively  shaUow,  average  slope  of  -0.034  for  the 
second  harmonic  and  -0.032  for  the  third  harmonic.  Again,  this  indicates  a  weak  dependence  of 
onset  distance  on  the  nonlinear  parameter,  a  medium  determined  property. 
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Figure  6-5  illustrates  the  effect  of  varying  path  attenuation  from  1x10  ^  to  10,  for  varying  B/A 
ratios,  on  onset  distance  of  second  and  third  harmonic  signals  from  a  fundamental  frequency 
signal  of  1.0  MHz,  10  ps  pulse  duration. 


Path  Attarualion  (d  BJb  mAIHz) 

Figure  6-5.  Effects  of  varying  path  attenuation  on  onset  distance  for  second  and  third  harmonic  signal  from  a  pulsed  sinusoid  of 
1.0  MHz,  10  ps  pulse  duration.  Note  that  as  path  attenuation  increases,  onset  distance  decreases  in  a  “pseudo-first  order 
manner  as  a  function  of  path  attenuation. 

As  the  results  in  Figure  6-5  suggest,  the  onset  of  the  third  harmonic  is  a  first  order  process, 
whereas  the  onset  distance  for  the  second  harmonic  appears  to  be  higher  order,  perhaps  second 
order.  Apparently,  there  are  at  least  two  competing  mechanisms  for  second  harmonic  onset 
distance.  The  first  mechanism  dominates  the  process  for  path  attenuation  less  than  5.00 
dB/cm/MHz  and  the  second  mechanism  dominating  for  path  attenuation  greater  than  5.00 
dB/cm/MHz.  Further  investigation  is  warranted  for  this  condition. 
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Effects  of  Varying  Input  Signal  Strength  on  Onset  Distance 

Figvire  6-6  shows  the  effects  of  varying  input  signal  strength  on  the  formation  of  harmonics  in 

simulated  biological  medium:  system  parameters  set  to  the  foUowing  values:  wave  velocity:  1500 
m/s;  signal  attenuation:  1x10’’  dB/cm/MHz;  B/A  ratio:  7.0;  signal  frequency:  1.75  MHz, 
observation  distance:  5.25  cm.  As  input  signal  strength  increases,  the  onset  distance  for  second 
and  third  harmonics  decreases.  The  dependence,  as  suggested  by  the  curvature  of  the  plots,  is 
most  likely  higher  than  first  order.  However,  if  first  order  is  assumed,  then  the  slopes  are 
indicated  by  the  dashed  lines:  -0.53 'and  -0.50  for  second  and  third  harmonics,  respectively. 


Figure  6-6.  Elfects  of  varying  input  signal  amplitude  on  harmonic  onset  distances.  Input  signal  strength  varies  from  0.5  MPa 
to  4.0  MPa.  System  parameter  settings:  wave  velocity:  1500  m/s;  frequency  dependent  attenuation:  1x10  dB/cm/MHz;  B/A 
ratio:  7.0;  signal  frequency:  1.75  MHz;  observation  distance:  5.25  cm. 

Effects  of  Varying  Input  Signal  Frequency  on  Onset  Distance 

Figure  6-7  shows  the  effects  of  varying  input  signal  frequency  on  the  formation  of  harmonics 


in  simulated  biological  medium:  system  parameters  set  to  the  following  values:  wave  velocity. 


1500  m/s;  signal  attenuation:  1x10’’  dB/cm/MHz;  B/A  ratio:  7.0;  signal  ampUtude:  0.5  MPa; 
observation  distance:  5.25  cm. 
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Figure  6-7  Effects  of  varying  input  signal  frequency  on  harmonic  onset  distance.  System  parameter  settings;  wave  velocity: 
1500  m/s;  frequency  dependent  attenuation:  Ixlff’  dB/cmMHz;  B/A  ratio:  7.0;  signal  input  amplitude:  0.5  MPa;  observation 
distance:  5.25  cm. 

The  plots  in  Figure  6-7  exhibit  slight  curvature  for  third  harmonic  onset  distance  dependence 


on  input  signal  frequency,  with  significant  curvature  for  second  harmonic  onset  distance 
dependence.  If  first  order  process  is  assumed,  then  the  slopes  are  indicated  with  the  dashed  lines; 
-0.45  and  -0.47  for  second  and  third  harmonics,  respectively. 


Narrowband  Sinusoid:  Image  Formation 

This  section  discusses  the  effects  of  parameter  variation  on  the  formation  of  ultrasoimd  B-scan 
images  using  the  same  computational  parameters  as  in  onset  distance  computation,  except  where 
the  particular  parameter  of  interest  is  varied. 

Effects  of  Frequency  Variation  on  Image  Generation 

Figure  6-8  shows  the  NUPROP  phantom  generation  module  (see  also  Figure  4-3).  The 

phantom  consists  of  three  ellipses.  The  first  ellipse  (Object  1)  is  the  largest  structure  (white)  and 
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is  positioned  as  the  center  of  the  256x128  array,  no  rotational  orientation.  The  second  ellipse 
(Object  2)  is  the  smallest  structure  (black),  is  rotated  by  ^5  degrees  (positive  angles  measured 
counterclockwise)  and  positioned  8  pixels  right  of  center,  10  pbcels  above  center.  The  third 
ellipse  (Object  3)  is  the  mid-sized  structure  (gray),  centered  within  the  256x128  array,  no 
rotational  orientation. 


Settings  for  Generating  Ultrasound  Phantom  Used  in  Figure  6-5  Sequence 


i|  NUPROP:  Tissue  Generalion  Module 


RIe  Execute 


BBC 


Tissue  SiTniation 


Parameters  Backgiounch 


B  ackground  Variation  Level  1  •  e-1 


Paianneters  Object  1 

ParaffTieters  2 - 

Pararrieters  Object  3  . . 

NorninalLevelllO-  | 

NoiTiinaiLevel|‘9-5  | 

Nominal  Level|'3-5  j 

Variation  Level|2.e-1  | 

Variation  Lcvdl^-e-l  [ 

Variation  Levd]  3-6-1  j 

Major  AwsModifwl  1-26-1  | 

Major  A»sModifier|  5-6-2  | 

Major  Modifief|l  -26-1  | 

Minof  Axis  ModJief|3'5e-1  j 

Minor  Axis  Modiier|7.e-2  | 

Minor  Axis  Modfia|3.5e-1  j 

Object  Ewtenl  Modifo  1 9  e-T  | 

0  bject  EHteht  Modifier|  9-^*1  | 

Object  Extent  Modlffef|7  e-l  | 

Object  Rotalbn  Angle  |0  | 

Object  Rotation  Angte| '45.  | 

Object  Rotation  AnglejO  | 

□  bject  S  hif  1  in  X-Direclion|  0  | 

0  biect  Shn  in  X-D  irection|9  | 

□biertShftinX-Direct»n|0  | 

Object  St#t  in  Y-Diection|0  | 

Object  Shift  in  Y*Diredion|‘iO  [ 

Object  SMt  in  Y-Directwn|0  | 

Pristine  Simi^ation 
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Figure  6-8.  NUPROP  Tissue  Generation  Module.  Settings  used  to  generate  simulated  ultrasound  images  for  sequence  in 
Figure  6-9.  Three  elliptical  objects  were  generated,  forming  the  complex  tissue  object.  Object  1  is  the  outer  ellipse  (white). 
Object  2  is  the  small  ellipse  at  a  45  degree  angle  (black).  Object  3  is  the  inner  ellipse  (gray). 

The  tissue  phantom  depicts  an  ovoid  tissue  structure  possessing  a  thick  membranous  wall,  a 
tissue  filled  core,  and  a  small  cyst  in  relative  positions  shown  in  the  pristine  simulation.  This 
module,  with  the  settings  shown  in  Figure  6-8,  was  used  to  generate  the  ultrasound  phantom  that 
produced  the  image  sequence  in  Figure  6-9. 
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Figure  6-9.  Effects  of  varying  center  frequency  of  a  narrow-band  sinusoidal  ultrasound  pulse  on  simulated  image  formation. 
Conditions  that  generated  this  sequence  of  images  are  specified.  Frequency  was  varied  from  1.0  MHz  to  3.0  MHz. 
Fundamental  (F);  Second  Harmonic  (2);  Third  Harmonic  (3):  and  Phantom  (P). 

Figure  6-9  shows  that,  as  the  initial  center  frequency  (Fundamental)  increases  from  1.0  MHz 
to  4.0  MHz,  the  image  resolution  improves.  The  structure’s  membranous  wall  is  resolved  at  4.0 
MHz.  The  membranous  wall  is  also  resolved  at  the  2.0  MHz  second  harmonic  frequency.  This 
suggests,  for  this  particular  phantom,  that  the  membranous  wall  possesses  a  high  correlation  with 
4.0  MHz  signals. 

The  target  cyst  is  detectable  at  all  four  fimdamental  frequencies,  with  the  “best”  detectability 
at  3.0  MHz.  Detectability,  in  this  case,  is  defined  as  the  highest,  localized  intensity  “spot”  in  the 
image.  The  cyst  is  localized  in  the  harmonic  images  of  each  frequency  as  well,  with  the  “best” 
localization  occurring  at  1.0  MHz  second  harmonic  and  3.0  MHz  second  harmonic. 
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The  2.0  MHz  second  harmonic  image  (4.0  MHz)  shows  marked  improvement  in  longitudinal 
resolution  over  the  4.0  MHz  fundamental  image.  This  result  is  not  unexpected,  since  literature  is 
replete  with  actual  experimental  results  suggesting  this  improvement.  The  average  signal  strength 
of  the  second  harmonic  image,  from  1.0  MHz  to  4.0  MHz,  is  —14.3  dB.  The  average  signal 
strength  of  the  third  harmonic  image,  over  the  same  frequency  spread,  is  -21.4  dB. 

Effects  of  Varying  B/A  Ratio  on  Image  Generation. 

Figure  6-10  shows  the  effect  of  varying  B/A  ratio  on  image  generation  using  a  narrowband 

sinusoid  with  center  frequency  of  1.75  MHz,  initial  amplitude  of  0.5  MPa,  velocity  of  1500  m/s, 
and  frequency  dependent  attenuation  of  IxlO  ’  dB/cm/MHz.  The  phantom  used  is  shown  in  the 
upper  right  comer.  Images  associated  with  the  fundamental,  second  harmonic,  and  third  harmonic 
are  represented  in  columns  1,  2,  and  3,  respectively.  B/A  ratios  of  5,  6,  and  10  were  used  (rows 
1,  2,  and  3  respectively).  Images  associated  with  the  difference  between  B/A  ratio  5  and  6,  5  and 
10,  and  6  and  10  for  the  fundamental,  second  harmonic,  and  third  harmonic  simulations  are 
presented  in  rows  4,  5,  and  6,  respectively.  Difference  images  are  presented  in  dB  scale. 

The  B/A  ratio  has  very  little  effect  on  the  fundamental  image,  as  is  readily  apparent  in  the 
difference  images  in  Figure  6-10.  The  effect  is  most  prominent  for  the  third  harmonic  difference 
images,  with  the  second  harmonic  exhibiting  some  small  differences  as  well.  In  each  case  for  the 
fimdamental  and  second  harmonic  image  differences,  however,  the  differences  are  negligible. 
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Figure  6-10.  Effects  of  varying  B/A  ratio  on  simulated  ultrasound  images.  The  phantom  used  is  shown  in  upper  right  comer. 
Images  associated  with  fundamental,  second  harmonic,  and  third  harmonic  are  shown  in  columns  1,  2,  and  3,  respectively.  B/A 
ratios  of  5,  6,  and  10  are  rows  1,  2,  and  3,  respectively.  Difference  images  5  -  6,  5  -  10,  and  6  ~  10  are  shown  in  rows  4,  5, 
and  6,  respectively.  Difference  images  are  presented  in  dB. 


Figure  6-11.  Histogram  of  Difference  images  in  Figure  6-8.  The  difference  images  from  Figure  6-7  are  redisplayed.  The 
histogram  of  the  fundamental  image  differences,  second  harmonic  image  differences,  and  third  harmonic  image  differences  are 
presented.  The  third  harmonic  difference  image  histogram  shows  appreciable  distribution,  whereas  the  fundamental  and 
second  harmonic  difference  image  histograms  essentially  show  no  difference. 
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Figure  6-11  shows  the  histogram  of  the  difference  images,  fiirther  confirming  the  negligible 
differences  between  fundamental  images  and  second  harmonic  images,  with  appreciable 
differences  appearing  in  the  third  harmonic  difference  images.  It  should  be  noted  that,  for  third 
harmonic  images,  the  maximum  signal  levels  range  fi’om  -44.56  to  -40.77  dB  over  the  span  of 
B/A  ratios  presented  in  Figure  6-10. 

The  effect  of  varying  B/A  ratio  is  manifested  in  the  overall  dB  level  for  the  second  and  third 
harmonic  images.  As  the  B/A  ratio  increases,  the  relative  signal  strength  increases  as  shown  in 
Figure  6-12. 


Figure  6-12.  EflFect  of  varying  B/A  ratio  on  signal  strength  for  images  generated  using  1.75  MHz  center  frequency,  narrowband 
sinusoid  insonification:  signal  amplitude:  0.5  MPa;  observation  depth:  4.2  cm;  velocity:  1500  m/s;  no  appreciable  frequency 
dependent  attenuation. 

Figure  6-12  suggests  that  for  typical  biological  tissue  types,  the  second  and  third  harmonics 
will  require  at  least  -27.79  to  -23.92  dB  and  ^6.67  to  -39.79  dB  of  gain  compensation, 
respectively,  to  produce  comparable  ultrasound  images  as  shown  in  Figure  6-10,  given  the 
phantom  used. 
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Effects  of  Varying  Signal  Amplitude  on  Image  Formation 

Figure  6-13  shows  the  effects  of  varying  the  input  signal  amphtude  on  image  formation. 

Signal  amplitudes  of  1.0,  2.0,  and  3.0  MPa  with  wave  velocity  of  1500  m/s,  signal  attenuation  of 
IxlO’  dB/cm/MHz,  observation  distance  of  2.5  cm.  The  images  generated  at  each  frequency 
were  pair-wise  differenced:  1.0  —  2.0  MPa,  1.0  —  3.0  MPa,  and  2.0  —  3.0  MPa.  As  can  be  readily 
seen,  signal  amplitude  has  little  apparent  visual  effect  on  the  fundamental  and  second  harmonic 
difference  images,  as  is  confirmed  by  the  difference  images’  histogram  plots.  The  primary  effect  is 
in  the  overall  signal  strength,  in  dB.  The  second  harmonic  images  have  signal  strengths  of -19.6 
dB,  -22.0  dB,  and  -23.0  dB  for  1.0, 2.0,  and  3.0  MPa  input  signal  amplitudes,  respectively. 

Fundamental  Second  Harmonic  Third  Harmonic 

1.0 

2.0 

3.0 

-  2.0 

3.0 
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Figure  6-13.  Effects  of  varying  signal  amplitude  on  image  formation.  Input  signal  strengths  of  1.0  MPa,  2.0  MPa,  and  3.0 
MPa  were  used.  System  parameters  were  set  as  follows:  wave  velocity:  1500  m/s;  B/A  ratio:  7.0;  observation  distance:  2.5  cm; 
frequency  dependent  attenuation;  1x10'*  dB/cm/MHz. 
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Differences  in  images  can  be  detected  at  the  third  harmonic,  where  the  histogram  plot  reveals 
at  the  high  and  low  ends  of  the  intensity  levels,  as  evidenced  in  Figure  6-14.  Again,  the  primary 
difference  between  the  three  input  signal  strength  images  is  manifested  in  the  overall  signal 
strength:  -29.3  dB,  -39.5  dB,  and  -40.3  dB  for  1.0,  2.0,  and  3.0  MPa  input  signal  strengths, 
respectively. 

Difference  Images  for  Amplitudes  of  1 .0, 2.0,  and  3.0  MPa 
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Figure  6-14.  Difference  images  and  histograms  associated  with  third  harmonics  for  1.0,  2.0,  and  3.0  MPa  signal  amplitudes. 
Differences  between  fiindamentals  and  second  harmonics  were  negligible.  System  parameters  were  set  as  follows:  wave 
velocity:  1500  m/s;  B/A  ratio:  7.0;  observation  distance:  2.5  cm;  frequency  dependent  signal  attenuation:  1x10"^  dB/cm/MHz. 

The  histogram  in  Figure  6-14  suggests  that  the  differences  are  primarily  due  to  noise  at  the 
lower  end  of  the  pixel  value  distribution,  whereas  differences  at  the  high  end  of  the  pixel 
distribution  may  be  attributable  to  actual  image  differences.  In  either  case,  the  range  of  values 
where  little  to  no  differences  occur  in  the  third  harmonic  difference  images,  namely  15  to  230  in  a 
range  from  0  to  255,  suggests  that  these  differences  are  negligible  as  well. 
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Effects  of  Varying  Wave  Velocity  on  Image  Formation 

Figure  6-15  shows  the  effects  of  varying  wave  velocity  on  the  formation  of  fundamental, 

second,  and  third  harmonic  simulated  images  using  the  phantom  (cyst  image)  generated  by 
NUPROP  (see  Figure  6-8).  The  system  parameters  were  as  follows:  input  signal  amplitude:  0.5 
MPa;  B/A  ratio:  7.0,  observation  distance:  2.5  cm;  frequency  dependent  attenuation:  1x10 
dB/cm/MHz.  Wave  velocity  was  varied  from  1400  to  1600  m/s,  in  100  m/s  increments.  As  is 
apparent  from  the  images  in  Figure  6-15,  there  is  very  Uttle  visual  difference  between  fundamental 
images  (1400  -  1500  m/s,  1400  -  1600  m/s,  and  1500  -1600  m/s)  and  second  harmonic  images 
(fundamental  at  1400  minus  fundamental  at  1500  m/s,  fundamental  at  1400  minus  fundamental  at 
1600  m/s,  and  fundamental  at  1500  minus  fundamental  at  1600  m/s)  and  second  harmonic  images 
(second  harmonic  at  1400  minus  second  harmonic  at  1500  m/s,  second  harmonic  at  1400  minus 
second  harmonic  at  1600  m/s,  and  second  harmonic  at  1500  minus  second  harmonic  at  1600  m/s). 

Differences  are  apparent  at  the  third  harmonic  as  illustrated  in  Figure  6-16.  Although  the 
actual  difference  images  do  not  exhibit  much  variation,  the  histograms  of  the  difference  images 
reveal  significant  difference  in  third  harmonic  at  1400  m/s  minus  at  1500  m/s  third  harmonic,  third 
harmonic  at  1400  m/s  minus  third  harmonic  at  1500  m/s,  and  third  harmonic  at  1500  m/s  minus 
third  harmonic  at  1600  m/s.  The  differenees  occur  at  the  periphery  of  the  distributions, 
suggesting,  at  the  lower  end  of  pixel  values,  the  differences  are  approaching  the  noise  floor  of  the 
proeess,  whereas,  at  the  high  end  of  the  pixel  value  distribution,  the  variation  is  due  to  actual 
differences  in  image  eontent.  So,  over  a  relatively  broad  range  of  pixel  values,  approximately  40 
to  180  in  a  range  from  0  to  255,  the  differences  in  third  harmonic  images  for  the  combinations  of 


the  three  velocities  is  negligible. 


•53 

JO 

> 


a> 

o 

$=l 

a> 

M 

ta 


Fimdamerital  S  econd  Harm  onjc _ TliiidHaiin^]i£. 


1400 


1600 


1800 


1400  -  1600 


1400  -  1800 


1600  -  1800 


Figure  6-15.  Effects  of  varying  wave  velocity  on  image  formation.  As  apparent  in  the  images,  very  little  difference  exists 
between  the  fimdamental  and  second  harmonic  images  of  varying  velocities.  Differences  are  manifested  in  the  third  harmonic 
image. 


Difference  Images  for  Velocities  of  1400, 1500,  and  1600  m/s 


Figure  6-16.  Difference  images  and  histograms  associated  with  third  harmonics  for  1400,  1500,  and  1600  m/s  wave  velocities. 
Differences  between  fundamentals  and  second  harmonics  were  negligible. 
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The  Apodized  Chirp:  Experiment  versus  NUPROP  modeling 

Up  to  this  point,  discussions  have  either  been  concerned  with  single  frequency  sinusoidal  input 
signals  to  determine  fundamental  eharacteristics  such  as  harmonic  onset  distances  and  beam 
profiles  at  various  penetration  depths;  or  with  narrowband  (multi-frequency)  sinusoidal  input 

signals  typically  used  for  tissue  imaging  (B-Scan  formation). 

One  novel  approach  is  the  use  of  a  linear  frequency  modulated  sinusoid  or  the  chirp. 
Appendk  III  discusses  NUPROP’s  implementation  of  nonlinear  chirps  and  their  effieacy.  This 
section  compares  NUPROP  simulation  of  an  apodized  chirp  with  experimental  data. 

Figure  6-17  shows  the  basic  frequency  modulated  sinusoid  or  chirp  signal.  Each  signal 
represented  exhibits  the  characteristic  “sawtooth”  structure  indicative  of  harmonics  present  in  the 
waveform.  Note  also  that  the  “real”  signals  have  a  higher  maximum  positive  value  than  maximum 
negative  value  whereas  the  “imaginary”  signals  are  symmetric  about  zero  amplitude.  Even  so,  the 
mean  amplitude  value  for  either  “real”  or  “imaginary”  signals  is  zero. 

Each  harmonic  frequency  within  the  propagated  chirp  signal  adds  “in-phase”  such  that  the 
zero  crossings  remain  fixed  throughout  the  harmonic  generation  process.  Another  indicator  that 
the  harmonics  are  adding  “in-phase”  is  that,  within  each  chirp  cycle,  the  waveform  is  sawtoothed. 
In  other  words,  if  any  appreciable  dispersion  existed  within  the  propagated  chirp  signal,  the 
coherent  nature  of  the  waves  would  necessarily  possess  destructive  interference,  eausing  the 
waveform  to  be  severely  disrupted.  Figure  6-18  illustrates  this  point.  The  smaller  circled  areas 
show  where  the  zero  crossings  are  changed.  The  larger  circled  area  shows  where  the  signal 
exhibits  destructive  interference. 
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Figure  6-18.  Illustration  of  the  effects  of  dispersion  or  phase  shifting  and  destructive  interference  on  coherently  added 
harmonic  signals.  Dispersion  is  indicated  within  the  small  circles.  In  this  example,  a  small  phase  shift  was  added  to  each 
harmonic  signal  causing  the  harmonics  to  eventually  add  out  of  phase,  manifested  in  the  destructive  interference  at  the  higher 
frequencies,  depicted  in  the  larger  circle. 
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TOP  VIEW 


SIDE  VIEW 


Figure  6-19.  Experimental  apparatus  used  to  collect  chirp  signal  data.  The  apparatus  consists  of  a  2.5  MHz  transducer 
acoustically  coupled  to  the  wall  of  the  water  tank.  The  hydrophone  is  attached  to  a  X-Y-Z  electronically  actuated  stage  for 
precise  positioning  within  the  water  tank  (9.5  cm  on-axis  from  the  face  of  the  transducer).  The  signal  from  the  hydrophone  and 
the  signal  generator  are  displayed  on  an  oscilloscope.  The  signals  are  digitized  and  sent  to  a  computer  for  storage  (digitizer  and 
computer  not  shown) 

Figure  6-19  is  a  cartoon  rendition  of  the  apparatus  used  to  collect  the  one  way,  on  axis  wave 
propagation  data.  The  transducer  used  was  a  Type  280  2.5  MHz  circular  transducer  (1.5  cm 
diameter).  The  hydrophone  used  was  a  Medisonics  Mark  II,  gain  set  to  zero.  A  Polynomial 
Wave  Synthesizer  Model  2020  was  used  to  generate  the  input  chirp  waveform.  The  equation  to 
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where  Mis  1x10*  and  |j  is  IxlO"*. 


The  signal  from  the  hydrophone  was  displayed  using  a  HAME  13,  60  MHz  Model  HM604 
Oscilloscope:  x-settings:  0.1  ps  delay,  time  division:  10  ps;  y-settings:  20  mV/0.1  division.  Data 
was  digitized  using  a  Data  Precision  Model  6100B  with  the  following  settings. 


x:  offset:  5.40  ps;  scale:  xl 


y:  offset:  O.OQv;  scale:  x2;  ±  500  mV;  AC  coupled 
timebase:  256  to  400  points;  period:  27.7  ns;  -1.662  ps  delay 


NUPROP  simulation  settings  were  as  follows: 


Parameter 

Value 

Initial  Frequency, 

1.5  MHz 

Sweep  Rate,  sweep 

2.9x10'^ 

Wave  Velocity,  c© 

1500  m/s 

Nonlinear  Parameter,  0 

5.0 

Initial  Wave  Amplitude,  Uo 

3.5  MPa 

Rayleigh  Distribution  Factor,  cc 

1.1x10* 

The  shock  parameter  distance,  I,  was  calculated  using  equation  3-8.  Each  distance  dependent 


amplitude  modifier  was  calculated  using 
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2’U^-L 

ri'Z 


r 


ri'Z 


(6-4) 


where  n  is  the  frequency  subscript  (i.e.  fimdamental  equals  1,  second  harmonic  equals  2,  etc.)  and 
z  is  the  observation  distance. 


Each  constituent  chirp  was  computed  using 


U.=A,-U,-e 


-i2nn{fo+sweep‘tyt 


(6-5) 


The  experimental  data,  waveform  at  the  observation  distance  of  9.5  cm,  exhibited  a  skewed 
distribution  introduced  by  the  transducer.  So,  an  apodization  function  based  on  the  Rayleigh 


distribution  was  used  such  that 

Apodization  =  STEP(t) 

a 


(6-6) 


where  STEP(t)  is  defined  as 


STEP(t-t,)  =  \ 


1,  for  t>tg 
0.5,  for  t  =  t„ 
0, for t  <t„ 


(6-7) 


with  to  equal  to  zero  for  this  example.  The  final  relationship  for  each  chirp  constituent  then  was 
computed  using 


CHIRP.  =  Apodization -U^ 


(6-8) 
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To  generate  the  composite  signal,  each  constituent  chirp  signal,  U„,  was  coherently  added. 


producing  the  results  in  Figure  6-20.  The  near  overlay  of  NUPROP’s  simulation  with 
experimental  data  confirms  the  hypothesis  that  each  fi-equency  within  a  given  chirp  sweep  adds 
coherently  and  “in  phase”  to  produce  the  sawtoothed  waveform  indicative  of  the  presence  of 
harmonics.  AdditionaUy,  the  hypothesis  that  the  chirp  could  be  substituted  as  the  bases  fimctions 
for  the  Fubini  solution  using  the  same  amplitude  weighting  Bessel  functions  is  confirmed. 
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Time  Domain  Representation  Fourier  Domain  Representation 


Solid  Lines:  Experimental  data 
Gray  Lines:  NUPROP  Simulation 

Figure  6-20  Comparison  of  NUPROP’s  simulated,  apodized  chirp  signal  and  experimental  data.  Gray  trace  is  NUPROP 
simulation,  black  trace  is  experimental  data. 

Figure  6-21  shows  the  reconstructed  constituents  of  the  chirp  signal  from  experimental  data. 
Each  constituent  within  the  experimental  data  was  reconstructed  from  correlations  between  the 
experimental  data  (penetration  depth  of  9.5  cm)  and  NUPROP’s  apodized  constituents  (1.5  MHz 
initial  frequency  for  the  fundamental,  3.0  MHz  initial  frequency  for  the  second  harmonic,  and  4.5 
MHz  initial  frequency  for  the  third  harmonic)  at  a  penetration  depth  of  9.5  cm.  As  is  evident  in 
Figure  6-21,  the  “fundamental”  chirp  (initial  frequency  1.5  MHz)  dominates  the  signature,  with 
sufficient  contributions  from  higher  harmonics  to  produce  the  “sawtoothed”  nature  of  nonlinearly 
propagated  waves. 
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Figure  6-21.  Reconstruction  of  experimental  data’s  chirp  constituents.  Each  constituent  chirp  is  produced  by  coirelating  the 
experimental  data  with  a  NUPROP  generated  apodized  chirp:  fundamental  chirp  initial  frequency,  1.5  MHz  with  a  sweep 
rate  of  2.7xl0'^,  second  harmonic  chirp  beginning  at  twice  and  the  third  harmonic  at  three  times  To. 

Simulation:  Narrowband  Sinusoid  versus  Apodized  Chirp 

This  section  presents  a  comparison  between  the  narrowband  sinusoidal  input  signal  and  the 
apodized  chirp  for  image  reconstruction,  simulations  only.  Figure  6-22  compares  a  narrowband 
sinusoid  and  an  apodized  chirp  signal.  The  center  frequency  of  the  narrowband  sinusoid  is  chosen 
to  match  the  average  frequency  of  the  chirp.  Additionally,  the  apodization  of  the  narrowband 
sinusoid  and  the  chirp  produce  matched  bandwidths. 
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Narrowband  Sinusoid 


Chirp 


1:  in  WDOndS 


Fundamental 

2nd  harmonic 

3rd  harmonic 

Phantom 
S/N:  5x10^ 


Figure  6-22.  Comparison  of  narrowband  sinusoid  and  apodized  chirp  signal  simulate  B-scan  images  with  signal  to  noise  ratio 
ofSxlO’. 

There  is  a  slight  gain  improvement  when  comparing  the  narrowband  sinusoid  second  harmonic 
image  to  its  chirp  counterpart,  approximately  2  dB  gain  improvement.  Though  this  may  seem 
significant,  the  reconstructed  B-Scan  images  reveal  that  the  narrowband  sinusoid  produces  a 
better  quality  image,  where  quality  is  measured  by  two  properties:  1)  preservation  of  structural 
integrity,  where  the  harmonic  image  shows  strong  indicators  of  structure  present  in  the  phantom, 
and  2)  detection  of  significantly  dijfferent  object  types  within  the  phantom,  object  discrimination. 
The  second  harmonic  image  for  the  narrowband  sinusoid  input  exhibits  strong  signature  returns 
fi-om  the  cyst  object  within  the  phantom.  The  second  harmonic  image  for  the  apodized  chirp  does 
not  show  such  indicators.  In  fact,  the  chirp  images  tend  to  look  “blurred,”  without  strong 
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structural  integrity,  nor  little  to  no  detection  between  the  larger  ovular  shapes  and  the  cyst.  So, 
for  relatively  high  signal  to  noise  ratios,  the  narrowband  sinusoid  tends  to  produce  “better  quahty” 

B-Scan  images. 

The  question  arises,  “when  does  the  chirp  produce  a  reasonably  significant  gain  in  image 
quality  vcrsus  the  naiTOwband  sinusoid?”  Signal  to  noise  ratio  plays  a  crucial  part  in  answering 
this  question.  A  block  diagram  of  signal  noise  sources  is  illustrated  in  Figure  6-23. 


Figure  6-23.  Block  diagram  of  potential  additive  noise  sources  for  an  ultrasound  imaging  system.  Noise  can  be  intr^uced  at 
each  electronic  device  within  the  image  reconstruction  process:  transmit  device,  sensing  device,  oscilloscope,  and  digiti2er. 

NUPROP  models  noise  as  a  single  additive  quantity  given  by 

^ total  ~  ^xmit  '^^o-scope  digitizer 

The  noise  fi’om  the  oscilloscope  is  minimal  whereas  the  transmit,  receive,  and  digitizer  noise  are 
significant. 

Figure  6-24  illustrates  the  effect  of  decreasing  signal  to  noise  ratio  on  the  reconstruction  of 
images  using  narrowband  sinusoidal  and  apodized  chirp  input.  Signal  to  noise  ratio  ranges  fi’om 
5x10’  to  750. 
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Phantom 


Figure  6-24.  Effects  of  signal  to  noise  ratio  on  image  formation  for  narrowband  sinusoid  and  apodized  chirp.  Noise  source: 
signal  electronics  at  the  input  and  output. 

As  noise  increases,  the  ability  for  the  narrowband  sinusoid  to  form  a  B-Scan  image  decreases 
whereas  the  apodized  chirp  continues  to  form  an  image  throughout  the  signal  to  noise  ratio  range 
investigated.  These  results  suggest  that  the  apodized  chirp  can  potentially  image  structures  under 
very  poor  imaging  conditions:  objects  deeply  buried  or  highly  occluded  objects. 

Figure  6-25  shows  the  effect  of  increasing  the  chirp  sweep  rate  on  B-Scan  image 
reconstruction.  Increasing  the  sweep  rate  while  maintaining  the  apodization  constant  effectively 
increases  the  time-bandwidth  produet.  The  center  frequency  of  the  narrowband  sinusoid  had  to 
be  increased  to  match  the  average  frequency  of  the  apodized  chirp  to  make  the  comparison. 
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Narrowband  Sinusoid 


Chirp  w/  increased  sweep  rate 


Fundamental 

2nd  harmonic 

3rd  harmonic 

Phantom 
S/N:  5x10^ 


Figure  6-25.  Effect  of  increasing  chirp  sweep  rate  while  maintaining  appropriate  apodization  and  comparable  center  to  average 
frequencies. 

The  increase  in  sweep  rate  produces  improved  apodized  chirp  B-Scan  images  (fundamental, 
second  and  third  harmonics),  with  structural  integrity  and  discrimination.  These  results  suggest 
that  the  apodized  chirp  possesses  additional  flexibility  in  design  over  the  narrowband  sinusoid. 


Chapter  6:  Summary 

The  concept  of  onset  distance  is  introduced.  Onset  distance  is  the  propagation  distance 
associated  with  the  first  time  a  frequency  component  reaches  five  percent  of  its  maximum  value. 
From  the  computations  presented,  the  formation  of  harmonics  is  not  strongly  dependent  on  B/A 
ratio.  Figures  6-3  and  6-4  tend  to  support  Christopher’s  contention  that  the  B/A  ratio,  >9,  is  not 
significant  to  nonlinear  wave  formation  results,  with  slopes  of  -0.034  and  -0.032  for  second  and 
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third  harmonics,  respectively.  Figure  6-5  suggests  a  stronger  dependence  on  path  attenuation 
than  on  B/A  ratio,  with  slopes  of  -0.1030  and  -0.0916  for  second  and  third  harmonics, 

respectively. 

Figure  6-6  suggests  that  onset  distance  is  strongly  dependent  on  input  signal  amplitude,  where 
the  second  and  third  harmonic  onset  distance  “rate  constants”  are  -0.53  and  -0.52,  respectively. 
Likewise,  Figure  6-7  shows  a  similar  dependence  on  signal  frequency  is  given  as  -0.45  for  the 
second  onset  distance  dependence  and  -0.47  for  the  third  harmonic  onset  distance  dependence. 

When  varying  input  frequency,  image  formation  is  significantly  different,  as  evidenced  in 
Figure  6-9.  As  the  input  frequency  increases  from  1  MHz  to  4  MHz,  different  parts  of  the 
simulation  phantom  are  accentuated.  For  instance,  the  “structural  wall”  is  highlighted  at  the 
second  harmonic  of  2  MHz  and  at  the  fundamental  of  4  MHz,  suggesting  that,  for  this  particular 
phantom,  the  “structural  wall”  contains  a  fair  amount  of  “4  MHz”  energy  to  return  to  the  sensing 
transducer.  Also  apparent  is  the  “improved”  image  resolution  at  all  second  harmomc  frequencies, 
improved  over  the  fundamental  of  each  frequency  tested.  The  mottled  image  quality  is  typical  of 
coherent  imaging  system  output. 

When  varying  B/A  ratio,  input  signal  amplitude,  or  wave  velocity,  there  seems  to  be  very  little 
effect  on  image  formation  at  the  fundamental  and  second  harmonic  frequencies  as  evidenced  in 
Figures  6-10,  6-11,  6-13,  6-14,  6-15,  and  6-16.  Differences  appear  to  become  significant  at  the 
third  harmonic  level.  The  primary  effect  of  varying  these  parameters  is  manifested  in  the  overall 
signal  level  for  each  harmonic.  Figure  6-12  shows  that  the  signal  level  increases  as  the  B/A  ratio 
increases,  suggesting  a  “rule  of  thumb”  for  setting  gain  levels  for  various  material  types  to  be 
tested. 
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Results  in  Figure  6-20  shows  that  NUPROP’s  chirp  model  closely  approximates  the 
experimental  data  collected  using  the  apparatus  in  Figure  6-19,  using  a  Rayleigh  distribution 
apodization  function.  The  results  also  indicate  that  the  original  hypothesis  that  the  chirp  functions 
can  be  substituted  as  the  bases  functions  in  the  Fubini  solution  to  nonlinear  wave  propagation. 

The  results  of  comparing  the  narrowband  sinusoid  and  apodized  chirp,  Figures  6-22  through 
Figure  6-25  suggest  that,  under  typical  ultrasound  imaging  conditions,  the  chirp  provides  very 
little  image  improvement.  Only  when  the  conditions  deteriorate  significantly  will  the  chirp  s 
efficacy  become  apparent:  detecting  structures  deeply  buried  or  highly  occluded. 
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Chapter  7.  Chapter  Summaries  and  Recommendations 

This  chapter  compiles  all  the  chapter  simimaries  and  provides  recommendations  for  future 
work  related  to  this  dissertation. 

Chapter  Summaries 

From  Chapter  2,  the  literature  is  replete  vwth  examples  of  ultrasound  modeling  techniques  that 
capture  aspects  of  the  overall  process.  For  nonlinear  wave  propagation,  two  models  in  literature 
are  prominent:  Christopher’s  NLP  and  the  Khokhlov-Zabolotskaya-Kuznetsov  (KZK) 

Implementation.  Christopher’s  NLP  is  based  on  solutions  to  Burgers’  equation  and  is 
implemented  n^ing  a  C-based  code.  Averkiuo’s  Implementation  is  based  on  the  KZK  equation, 
that  includes  diffraction,  attenuation,  and  nonlinearities  in  a  single,  third  order  partial  differential 
equation. 

For  linear  propagation,  two  methods  are  presented:  angular  spectrum  and  Lommel 
formulation.  Angular  spectrum  implements  a  planar  geometry  derived  mathematical  convention 
using  a  propagation  factor  that  accounts  for  complex  field  changes  over  small  incremental 
distances.  Angular  spectrum  suffers  fi:om  computational  size  limitations,  thus  reducing  its  angular 
resolution  capabilities.  Also,  since  angular  spectrum  is  developed  using  planar  calculus,  focused 
beam  geometries  are  not  easily  accommodated. 

Angular  spectrum  is  applicable  throughout  the  propagation  path,  fi'om  just  in  fi'ont  of  the 
transducer  to  as  far  as  the  user  wished  to  compute  the  field.  Additionally,  angular  spectrum  is 
suitable  for  computing  field  propagation  fi-om  non-radially  symmetric  geometries. 

Lommel  formulation  implements  planar  as  well  as  focused  beam  transducer  geometries. 
However,  Lommel  is  applicable  for  circularly  symmetric  geometries  only.  Additionally,  Lommel 
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is  strictly  a  linear  propagation  process.  Therefore,  any  process  that  introduces  nonlineanties,  i.e. 
generation  of  harmonics,  will  necessarily  be  inaccurately  predicted  by  Lommel. 

Ultrasound  image  generation  models  dependent  strongly  on  the  tissue  characterization  or 
tissue  construction  model  used.  One  model  has  found  considerable  fevor  in  literature;  Bamber 
and  Dickinson’s  continuous  array  of  point  scatterers.  If  the  tissue  eharacteristic  of  interest  can  be 
modeled  as  an  array  of  continuous  point  scatterers,  then  Bamber’s  techmque  is  reasonably 
accurate  in  predicting  linear  ultrasound  effeets.  Bamber’s  techmque  involves  the  eonvolution  of 
the  tissue  array  with  the  point  spread  funetion  of  the  ultrasound  pulse.  Bamber  aceonplishes  the 
convolution  by  taking  the  process  to  the  Fourier  domain  and  multiplying  the  Fourier  transforms  of 
the  tissue  array  and  the  point  spread  function,  then  inverse  Fourier  transforming  the  product. 

Though  models  exist  for  the  various  aspects  of  the  ultrasound  process,  no  single  model 
encompasses  all  aspects  of  ultrasound  propagation  and,  ultimately,  ultrasound  imaging.  A  new 
model.  Nonlinear  Ultrasound  Propagation  (NUPROP),  is  introduced  that  captures  aspects  of 
nonlinear  wave  generation  and  propagation,  linear  wave  generation  and  propagation,  and  image 
generation. 

From  Chapter  3,  NUPROP  is  a  fairly  complex  simulation  and  analysis  tool  for  ultrasoimd 
modeling.  NUPROP  represents  an  end-to-end  modeling  environment  for  studying  ultrasound 
wave  propagation  and,  ultimately,  B-Scan  image  formation.  This  chapter  diseussed  the 
generation  of  objects  called  phantoms  using  Bamber’s  continuous  point  scatterer  concept.  The 
tissue  characteristic  of  interest  was  acoustic  impedance.  Input  signal  waveforms  used  in 
NUPROP  are  either  sinusoidal  or  frequency  modulated  sinusoidal,  either  apodized  (for  image 
formation)  or  non-apodized  (single  frequeney  or  continuous  used  for  code  verification  and 
harmonic  signal  analyses). 
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Nonlinear  propagation  was  modeled  using  either  the  Burgers’  equation  or  Fubini  solution 
approaches.  Blackstock  developed  a  closed  form  solution  to  the  Burgers’  equation  that 
possessed  the  shock  distance  parameter,  L,  which  is  calculated  from  system  parameters  of 
velocity,  initial  wave  amplitude,  nonlinear  parameter  {fi),  and  fundamental  frequency.  Likewise, 
the  Fubini  approach  possessed  this  same  shock  distance  limitation. 

Linear  propagation  was  modeled  using  Angular  Spectrum  and  Lommel.  Angular  spectrum 
has  the  advantages  of  being  applicable  throughout  the  propagation  path,  accommodate  non- 
radially  symmetric  transducer  geometries,  and  accurately  computes  nonlinear  or  harmonic 
generation.  Angular  spectrum  requires  a  trade-off  between  computational  speed/volume  and 
accuracy/precision  in  computing  trans-axial  beam  profiles.  Lommel,  on  the  other  hand,  computes 
the  on-axis  longitudinal  beam  profile  and  the  trans-axial  beam  profiles  at  depth  with  accuracy  and 
precision.  Unfortunately,  Lommel  formulation  suffers  from  two  major  limitations:  radial 
symmetry  only  and  linear  propagation  only.  Lommel  does  not  predict  harmonic  signal  formation 
nor  does  it  accurately  compute  off-axis  beam  profiles  for  non-radially  S5nnmetric  geometries. 
Lommel  does  not  accurately  compute  the  pressure  field  just  in  front  of  the  transducer.  A 
minimum  propagation  distance  is  required  before  Lommel  can  be  used. 

Image  generation  is  accomplished  using  the  Hilbert  transform,  which  converts  a  RF  signal  to 
an  “envelope”  detected  intensity  signal. 

Finally,  NUPROP  performs  some  rudimentary  post  processing  such  as  statistical  analyses, 
histogram  plotting,  and  simple  image  thresholding  based  on  the  statistics  of  the  images. 

From  Chapter  4,  NUPROP  is  a  modularly  design  IDL  code  that  models,  end-to-end,  nonlinear 
ultrasound  wave  propagation  and  ultrasound  image  formation.  NUPROP  is  based  on  Svidgets, 
or  controls,  in  the  terminology  of  some  development  environments.  Widgets  are  simple  graphical 
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objects  such  as  pushbuttons  or  sUders  that  aUow  user  interaction  via  a  pointing  device  (usuaUy  a 
mouse)  and  a  keyboard.  IDL  graphical  user  interfaces  are  constructed  by  combining  widgets  m  a 
tree-like  hierarchy.  Each  widget  has  one  parent  widget  and  zero  or  more  child  widgets.  There  is 
one  exception:  the  topmost  widget  (caUed  a  top-level  base)  is  always  a  base  widget  and  has  no 
parent.  So,  NUPROP  is  a  graphical  user  interface  (GUI)  providing  the  interactive  platform  for 
the  user  to  exercise  the  physics  associated  with  nonlinear  and  linear  wave  propagation. 

Menu  driven,  NUPROP  is  built  with  modules  (operations  that  are,  essentially,  stand-alone  m 
that  they  do  not  need  the  GUI  to  work)  and  functions  (modular  code  that  requires  parameter 
values  passed  through  the  GUI  to  operate).  The  complexity  of  each  module  or  function 
determines  the  computational  requirements,  hence  the  computational  time,  accuracy,  and 
precision. 

From  Chapter  5,  NUPROP’s  nonlinear  wave  generation  and  propagation  results  compared 
favorably  with  data  found  in  literature,  having  less  than  one  percent  difference  comparmg 
NUPROP’s  Fubini  implementation  with  existing  data  (Ryan,  [1963]  and  less  than  a  few  percent 
between  NUPROP’s  Burgers’  solution  and  existing  data.  NUPROP  also  accurately  predicted  on- 
axis  wave  amplitude  for  propagation  distances  comparable  to  Averkiuo’s  simulation  conditions. 
Additionally,  NUPROP  confirmed  Averkiuo’s  supposition  that,  for  equivalent  areas,  circular 
transducers  produce  equivalent  on-axis  longitudinal  amplitudes  as  rectangular  arrays  (P3-2 
geometry).  However,  NUPROP  also  showed  that  for  wave  or  field  amplitudes  off-axis, 
Averkiuo’s  assumption  is  invalid. 

In  comparing  trans-axial  beam  profiles  derived  fi*om  diffraction  coupled  with  nonlinear 
propagation,  NUPROP  compared  very  favorably  with  hterature  (Christopher,  [1997]),  having  less 
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than  a  few  percent  difference  in  trans-axial  beam  profile  amplitudes  for  2  MHz  and  its  second 
harmonic  at  4  MHz,  within  1 .5  cm  radial  distance  fi-om  the  axis. 

Comparing  NUPROP’s  RF  and  envelope  detection  simulated  images  with  Bamber’s  technique 
for  a  point  source.  Bamber’s  simulation  used  a  “skewed”  or  non-symmetric  temporal  pulse  shape, 
whereas  NUPROP  assumes  a  temporally  symmetric  beam  apodization.  Comparing  the  3  dB  point 
for  the  most  gradual  fall-off,  NUPROP  is  within  one  percent  in  computing  the  half  power  point, 
or  3  dB  down  fi'om  maximum  intensity. 

Given  the  favorable  comparison  of  NUPROP’s  simulations  of  nonlinear  wave  generation  and 
propagation  using  Fubini  or  Burgers’  partial  differential  equation  solutions,  linear  wave 
propagation  (diffraction)  using  angular  spectrum  or  Lommel  formulation,  and  image  formation 
using  Bamber’s  technique,  verification  of  NUPROP’s  viability  as  a  simulation  code  is  complete. 

From  Chapter  6,  the  concept  of  onset  distance  is  introduced.  Onset  distance  is  the 
propagation  distance  associated  with  the  first  time  a  fi*equency  component  reaches  five  percent  of 
its  mavimiim  value.  From  the  computations  presented,  the  formation  of  harmonics  is  not  strongly 
dependent  on  B/A  ratio.  Figures  6-3  and  6-4  tend  to  support  Christopher’s  contention  that  the 
B/A  ratio,  p,  is  not  significant  to  nonlinear  wave  formation  results,  with  slopes  of -0.034  and 
-0.032  for  second  and  third  harmonics,  respectively.  Figure  6-5  suggests  a  stronger  dependence 
on  path  attenuation  than  on  B/A  ratio,  with  slopes  of  -0.1030  and  -0.0916  for  second  and  third 
harmonics,  respectively. 

Figure  6-6  suggests  that  onset  distance  is  strongly  dependent  on  input  signal  amplitude,  where 
the  second  and  third  harmonic  onset  distance  “rate  constants”  are  -0.53  and  -0.52,  respectively. 
Likewise,  Figure  6-7  shows  a  similar  dependence  on  signal  frequency  is  given  as  -0.45  for  the 
second  onset  distance  dependence  and  —0.47  fi)r  the  third  harmonic  onset  distance  dependence. 
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When  varying  input  frequency,  image  formation  is  significantly  different,  as  evidenced  in 
Figure  6-9.  As  the  input  frequency  increases  from  1  MHz  to  4  MHz,  different  parts  of  the 
simulation  phantom  are  accentuated.  For  instance,  the  “structural  wall”  is  highlighted  at  the 
second  harmonic  of  2  MHz  and  at  the  fimdamental  of  4  MHz,  suggesting  that,  for  this  particular 
phantom,  the  “structural  waU”  contains  a  feir  amount  of  “4  MHz”  energy  to  return  to  the  sensing 
transducer.  Also  apparent  is  the  “improved”  image  resolution  at  aD  second  harmonic  frequencies, 
improved  over  the  fimdamental  of  each  frequency  tested.  The  mottled  image  quahty  is  typical  of 
coherent  imaging  system  output. 

When  varying  B/A  ratio,  input  signal  amplitude,  or  wave  velocity,  there  seems  to  be  very  little 
effect  on  image  formation  at  the  fimdamental  and  second  harmonic  frequencies  as  evidenced  in 
Figures  6-10,  6-11,  6-13,  6-14,  6-15,  and  6-16.  Differences  appear  to  become  significant  at  the 
third  harmonic  level.  The  primary  effect  of  varying  these  parameters  is  manifested  in  the  overall 
signal  level  for  each  harmonic.  Figure  6-12  shows  that  the  signal  level  increases  as  the  B/A  ratio 
increases,  suggesting  a  “rule  of  thumb”  for  setting  gain  levels  for  various  material  types  to  be 
tested. 

Results  in  Figure  6-20  shows  that  NUPROP’s  chirp  model  closely  approximates  the 
experimental  data  collected  using  the  apparatus  in  Figure  6-19,  using  a  Rayleigh  distribution 
apodization  function.  The  results  also  indicate  that  the  original  hypothesis  that  the  chirp  functions 
can  be  substituted  as  the  bases  functions  in  the  Fubini  solution  to  nonlinear  wave  propagation  is 
valid. 

The  results  of  comparing  the  narrowband  sinusoid  and  apodized  chirp.  Figures  6-22  through 
Figure  6-25  suggest  that,  under  typical  ultrasound  imaging  conditions,  the  chirp  provides  very 
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little  image  improvement.  Only  when  the  conditions  deteriorate  significantly  will  the  chirp  s 
efBcacy  become  apparent:  detecting  structures  deeply  buried  or  highly  occluded. 

Recommendations 

This  section  discusses  additional  experiments  that  could  be  conducted  to  further  vaUdate  the 
IDL  code  as  well  as  code  related  activities  and  potential  improvements  to  NUPROP.  Though 
couched  in  terms  of  NUPROP  code  validation  and  improvements,  these  recommendations  also 
serve  as  potentially  fiiiitful  research  areas  in  their  own  right  and  should  be  pursued  whether  or  not 
they  are  part  of  NUPROP  improvements.  Code  improvements  include  the  addition  of  the 
diffraction  integral  and  shift-and-add  methods  for  computing  linearly  propagated  wavefi-onts, 
complementing  angular  spectrum  and  Lommel;  addition  of  indices  commonly  used  in  climcal 
diagnoses,  namely  the  mechanical  index  which  is  an  indicator  of  cavitation  and  the  two  thermal 
indices  (cranial  and  soft  tissue),  which  are  indicators  of  thermal  stress  and  potentially  dangerous 
thermal  loading.  The  use  of  more  complex  phantoms  such  as  the  piecewise  homogeneous  model 
or  the  addition  of  phase  aberration  are  warranted. 

NUPROP  Code  Validation 

This  section  describes  potential  experimentation  to  vahdate  NUPROP  code  for  single 
fi-equency  sinusoid  cross-track  beam  profiling  and  pulse  sinusoid  image  formation  (B-Scan). 

Single  Frequency  Sinusoid 

NUPROP  represents  the  culmination  of  a  considerable  body  of  research  on  nonlinear  and 
linear  propagation  modeling.  Considering  nonlinear  propagation  first,  the  Fubini  solution  code 
has  been  verified  against  existing  data  and  results  fi-om  other  simulation  codes  and  validated 
against  experimental  data.  The  Burgers’  computational  model  has  been  verified  but  not  validated. 
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For  linear  propagation,  the  Lommel  formulation  module  has  been  verified  against  existing  data 
but  is  rudimentary,  at  best.  This  module  requires  more  work  to  include  angular  or  cross-track 
beam  formation  similar  to  the  angular  spectrum  module.  The  angular  spectrum  module  has  been 
verified,  but  not  validated  with  experimental  data.  Experiments  similar  to  the  apodized  chirp 
should  be  conducted  to  determine  the  off  axis  (cross-track)  beam  profiles  at  the  observation 
distance. 

Figure  7-1  illustrates  a  cartoon  of  the  experimental  apparatus  required  to  conduct  beam  profile 
determination. 


Figure  7-1.  Cartoon  rendition  of  beam  profile  measurement  apparatus.  Waveform  fi*om  the  signal  generator  drives  the 
transducer  (single  fi*equency  sinusoid).  The  hydrophone  is  accurately  positioned  using  the  computer  controlled  X-Y-Z  stage. 
The  Y-Position  is  fixed  as  the  X-Position  scans  in  micron  steps  across  the  ultrasound  beam.  Signals  fi-om  the  hydrophone  are 
displayed  on  an  oscilloscope,  with  further  processing  and  data  capture  using  a  signal  digitizer  and  computer  (not  depicted  in 
this  cartoon). 

The  transducer  is  driven  by  a  signal  fi-om  the  signal  generator:  a  single  frequency  sinusoid. 
The  signal  generator  signal  is  sent  to  an  oscilloscope  for  display.  The  hydrophone  is  positioned 
accurately  using  the  eomputer  eontroUed  X-Y-Z  stage.  The  Y-Position  is  fixed  while  the  X- 
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position  scans  in  micron  steps.  The  signal  from  the  hydrophone  is  sent  to  the  oscilloscope  for 
display.  Both  the  original  signal  and  the  hydrophone  signal  are  sent  to  a  digitizer  for  data  capture 
and  then  to  the  computer  for  processing  and  storage.  The  average  signal  strength  at  each  X- 
Position  is  computed  and  plotted  as  a  function  of  cross-track  displacement.  The  computed  cross¬ 
track  beam  profile  is  then  compared  with  NUPROP’s  predicted  cross-track  beam  profile  for 
validation. 

Pulse  Sinusoid  Image  Formation  (B-Scan) 

Figure  7-2  depicts  a  cartoon  rendition  of  the  experimental  apparatus  used  for  pulse  sinusoid 
image  formation  (B-Scan). 


to  osciltoscope 


Figure  7-2.  Illustration  of  experimental  apparatus  used  for  ultrasound  B-scan  imaging.  A  transducer  is  attached  to  a  single 
stage  computer  controlled  actuator  for  precision  placement  along  the  scan  direction.  A  wire  phantom  containing  several  wires 
spatially  displaced  is  insonified.  The  reflected  ultrasound  energy  is  captured  by  the  transducer.  The  B-scan  image  is  generated 
from  A-line  scans  of  the  wire  phantom  at  various  points  along  the  scan  path. 

Energy  introduced  by  the  transducer  at  each  position  interacts  with  structure  along  the 
propagation  path.  Energy  reflected  from  the  structures  in  the  propagation  path  is  collected  over 
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the  dwell  time  of  the  transducer.  The  time  signal  is  commonly  referred  to  as  an  A-Line  scan.  A- 
Line  scans  comprise  the  B-scan  image  after  envelope  detection. 


Code  Improvements 

NUPROP  currently  models  planar  wavefronts.  A  significant  inprovement  would  be  the 
introduction  of  focused  wavefronts.  This  can  be  accomplished  using  an  angular-spatial  phase 
filter  and  employing  angular  spectrum.  However,  this  method  is  not  efficient  as  pointed  out  by 
Liu  and  Waag  [1997].  The  introduction  of  the  diffraction  integral  method  or  the  more  empirical 
shift  and  add  method,  as  mentioned  in  the  Liu  and  Waag  [1997]  article,  are  more  efficient  and 
would  complement  the  Angular  Spectrum  and  Lommel  techniques  currently  incorporated  in 
NUPROP. 


Diffraction  Integral  Method 

The  diffraction  integral  method  relies  on  the  fact  that  the  pressure  field  satisfies  the  Helmholtz 
equation  such  that,  at  any  arbitrary  point,  the  value  of  the  pressure  field  can  be  expressed  as  a 


surface  integral 


(7-1) 


where  5^  is  an  arbitrary  surface  enclosing  the  point  r ,  r^  is  a  point  on  So,  and  is  the  outward 


pointing  normal  vector  of  So  at  r^ .  The  Green’s  fimction,  Girjrf),  is  a  general  solution  to  the 
Helmholtz  equation  such  that 
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G^lr^=g{R)  +  X(r) 


(7-2) 


where  R=r-r^  is  the  distance  between  r  and  r,, ,  g(R)  is  the  free-space  outgoing  Green’s 
function  as  defined,  and  xO")  is  an  arbitrary  function  satisfying  the  Helmholtz  equation 


(V^  ■^k^)x(y)  =  0  iit  the  volume  surrounded  by  the  surface,  So.  So,  the  field  at  an  arbitrary 
point  in  space  can  be  computed  using  the  field  at  a  known  point  in  space,  namely  on  the  surface 
surroxmding  the  point,  given  two  conditions:  1)  the  surface  is  comprised  of  a  plane  z=0  and  the 
hemisphere  at  infinity  enclosing  the  upper  half  space,  and  2)  the  behavior  of  the  pressure  field  at 
infinity  is  specified  to  satisfy  the  Sommerfeld  radiation  condition  in  Goodman  [1968].  The 


pressure  field  at  p(x,y,z)  is  given  as 


p{x,y,z)  =  ^  \\p{x„,y„S!i)  —-ikj-^dx„dy„ 


p{x,y,z)  =  -—\\ 


1  e'*'' 


dz„  R 


where  R  =  ^{x-x„y  ^^y-yof  • 

Shift-and-Add  Method 

The  Shifl-and-Add  method  can  be  derived  fi-om  the  diffraction  integrals  by  an  inverse  Fourier 
transform  of  (7-3)  with  respect  to  time  followed  by  a  discretization  of  the  integrals  in  space. 
Usually,  the  distance  R  is  sufiBciently  large  compared  to  the  wavelength,  X,  such  that 

ikz  ft  e'*"  ('7-4) 

Pa,i^,y,z)  =  -—\]Po,{Xo^yoSi)-^dx^(fy„ 


omitting  the  1/R  term  fi-om  (7-3)  and  signifying  that  the  equation  applies  to  a  single  fi-equency 
component,  o).  The  inverse  Fourier  transform  of  (7-4)  becomes 


±Q0 

y,hevep(x,y,z,t)  is  the  temporal  pressure  variation  at  point  (x,y,z)  and  p  denotes  the  differentiation 
ofp  with  respect  to  time.  The  subsequent  discretization  with  respect  to  x<,  andyo  yields 

1  (7-6) 

p{x,  y,Z,t)=  -r—  Ayo  ^  X  S  P^^om  ’  Ton  t  -  R„„  /  c)  2 

^TIC  m  n 

with  R„„  =  ^|ix-x„„f  +iy-yonf  • 

Both  the  Diffiaction  Integral  and  the  Shift-and-Add  methods  accommodate  focused 
wavefronts,  and  are  computationaUy  more  efficient  than  Angular  Spectrum  in  this  respect.  Since 
Fubini  provides  a  linear  approximation  of  a  nonlinear  process,  namely,  harmomc  signal  generation, 
coupling  nonlinear  and  linear  propagation  using  either  the  diffiaction  integral  or  shift-and-add 

methods  is  relatively  simple. 

Mechanical  Effects:  Cavitation  and  mechanical  Index  (MI) 

The  problem  of  cavitation  has  long  been  of  interest  in  hydrodynamics,  where  it  is  primarily 

responsible  for  the  destruction  of  propellers,  turbines,  and  other  moving  parts.  A  fundamental 

problem  is  understanding  the  mechanism  for  cavitation,  conditions  that  promote  the  onset  of 

cavitation.  As  the  name  suggests,  cavitation  is  the  formation  of  “holes”  in  the  liqmd,  usually 

under  tensile  stress.  Depending  on  the  conditions,  the  “holes”  may  be  filled  by  dissolved  gases  in 

the  Uquid.  However,  for  ideally  “pure”  Uquids,  cavitation  indicates  a  ripping  apart  of  the  structure 

of  the  liquid.  One  potential  measure  of  cavitation  is  the  mecharacal  index  (MI).  In  principle,  the 

use  of  the  mechanical  index  gives  the  user  a  means  of  judging  the  safety  of  clinical  applications  of 

ultrasound.  Unfortunately,  experimental  evidence  for  mechanical  damage  within  soft  tissues  from 
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diagnostic  ultrasound  is  limited,  at  best.  With  the  possible  exception  of  lung  hemorrhage,  it  is 
difficult  to  relate  values  of  MI  to  specific  hazards. 


Mechanical  Index  is  calculated  fi-om  the  measured  peak  negative  pressure,  p~ ,  in  MPa  as 


(7-7) 


where /is  the  fi-equency,  in  MHz.  Peak  negative  pressure  is  calculated  by  reducing  the  acoustic 
pressure,  measured  in  water,  using  an  attenuation  coefficient  of  0.3  dB  cm*  MHz*.  This 
coefficient  of  attenuation  is  an  approximate  model  of  biological  tissue,  assuming  homogeneous, 
single  layer  medium.  MI  indicates  the  amplitude  of  the  ultrasonic  pulses  being  used  at  any  time. 
Increased  pulse  amplitudes  result  in  proportionally  higher  MI  values.  The  rationale  is  that 
increased  amplitudes  will  potentially  cause  cavitation,  or  mechanical  damage.  Hence,  MI 
becomes  the  measured  parameter  for  determining  potential  mechanical  damage  in  biological 
tissue.  The  formation  of  harmonics  signals  produces  the  potential  for  very  high,  localized  velocity 
differentials.  Compression  and  expansion  within  a  very  small  distance  can  cause  severe  localized 
distortion.  Figure  7-3  illustrates  this  condition.  As  the  wave  becomes  more  and  more 
“sawtoothed”  due  to  harmonic  signal  formation,  the  wavefi-ont  can  exhibit  near  instantaneous 
tnavimiim  positive  to  maximum  negative  velocity  within  a  few  microns,  introducing  sufficient 
energy  to  break  chemical  bonds. 
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Nonlinear  Wavefront 


Figure  7-3.  Potential  wave  differential.  Region  within  a  few  microns  has  the  potential  of  experiencing  a  nearly  instantaneous 
velocity  differential  of  2.X . 

Nonlinear  effects  influence  the  achievable  peak  rarefactional  pressure  used  in  calculating  MI. 
This  peak  negative  pressure  tends  to  saturate  in  water  at  the  highest  diagnostic  outputs.  For 
unfocused  transducers,  Carstensen’s  [1996]  nonlinear  propagation  models  underestimate  the  peak 
negative  pressure  by  as  much  as  20%,  even  at  relatively  low  source  intensities.  Christopher’s 
nonlinear  model  [1997]  showed  an  underestimation  of  8%  for  the  same  transducer  studied. 
Christopher  employed  a  more  sophisticated  nonlinear  model  based  on  Burgers  equation.  The 
accuracy  and  executability  achieved  with  these  computational  techmques  are  critical  to  precise 
estimates  of  MI  in  vivo. 

The  source  intensity  tends  to  shift  the  beam  focus  closer  to  the  source,  even  in  “lossless” 
medium  of  water.  The  current  protocol  for  derating  climcal  applications  of  ultrasound  tends  to 
shift  the  maxima  closer  to  the  source.  Nonlinear  propagation  produces  a  strange  absorption 
coelBBcient  that  increases  in  magnitude  with  distance  fi-om  the  source.  As  a  consequence,  the  shift 
of  the  maximum  pulse  intensity  becomes  significant,  resulting  in  a  discontinuity,  or  potentially 


99 


damaging  mechanical  condition  within  the  tissue.  Even  with  this  potentially  damaging  condition, 
it  is  still  impossible  to  draw  general  conclusions  on  the  MI  error  magnitudes  introduced  by 
nonlinearities.  The  field  regime  in  which  these  discontinuities  exist  have  little  to  no  nonlinear 
effects.  Underestimates  as  great  as  a  factor  of  two  have  been  demonstrated  in  specific 
experimental  configmations. 

Thermal  Index  (TI) 

The  second  biophysical  mechanism  of  interest  is  tissue  heating,  “quantified”  using  thermal 
index  (TI),  an  estimate  of  the  tissue  temperature  rise  in  "C  which  might  occur  under  “reasonable 
worst-case  conditions.”  Three  thermal  indices  have  been  defined  in  literature:  1)  soft  tissue  index 
(TIS),  2)  bone  index  close  to  the  transducer  (Cranial  conditions,  TIC),  and  3)  bone  index  at  or 
near  the  transducer  focus  (TIB).  Thermal  index  is  defined  as  the  ratio  of  the  acoustic  power 
emitted  to  the  power  required  to  heat  a  particular  target  tissue  by  1  °C,  allowing  time  for  thermal 
equilibrium. 

Nonlinearities  will  limit  the  amplitude  of  fields  distant  fi-om  the  source,  but  they  also  increase 
local  absorption.  These  factors  are  oppose  each  other’s  biological  effects  giving  rise  to  nonlinear 
propagation  either  increasing  or  decreasing  localized  heating,  depending  on  which  mechanism  is 
dominant  under  the  observation  conditions. 

Thermal  indices  (TIs)  are  designed  to  estimate  spatial  maximum,  steady-state  temperature 
increments.  For  homogeneous  tissue,  the  thermal  index  is  defined  as  TIS  and  is  given  as  the 
temperature  increment  at  a  fixed  distance  of  ~2  cm  fi'om  the  source,  not  at  the  focus.  At  this 
distance,  nonlinear  propagation  enhances  local  thermal  increment  by  10%.  TIS  for  a  given 
transducer  configuration  is  expressed  by  one  of  two  mathematical  relationships 


100 


775  =  ^0  3 

210 


(7-8) 


where  hM^)  is  the  local,  derated,  temporal  average  intensity  in  mW/cm^  f  is  frequency,  or 


_  W,,{z)f 

210 


(7-9) 


where  WoM^)  is  the  total,  derated  acoustic  power  at  or  near  the  focus  where  the  effective  area  of 
the  beam  is  less  that  1  cm^. 

When  bone  is  the  target  tissue,  nonlinear  propagation  has  no  effect  on  its  absorption.  The 
only  effect  of  importance  is  the  reduction  of  local  fields  and,  thus,  reduction  of  thermal  loading  or 
heating.  Thermal  Index  for  bone  (TIB)  is  computed  using 


TIB  = 


50 


(7-10) 


where  z,  in  this  case,  is  the  depth  that  maximizes  TIB. 

Thermal  indices  used  in  clinical  applications  depend  on  total  acoustic  power,  tissue  absorption 
coefficient,  frequency,  beam  size,  tissue  thermal  capacity,  tissue  thermal  conductivity,  and 
estimates  of  blood  perfusion  for  tissue  cooling.  Thermal  index  computational  constants  (210  for 
TIS,  50  for  TIB)  are  derived  from  these  properties.  Since  the  propagation  distance  is  very  small 
for  determining  cranial  thermal  index  (TIC),  nonlinear  propagation  is  a  non-player. 

By  in  large,  nonlinear  wave  propagation  plays  a  minor  role  in  the  TI  Output  Display  Standard 
(ODS)  because  of  the  way  thermal  indices  are  defined  and  because  of  the  protocols  recommended 
for  measuring  acoustic  power.  In  contrast  to  mechanical  index  (MI),  thermal  indices  (TIs)  were 
developed  using  well-understood  heat  generation  and  difiusion  theories.  The  theories  are 
straightforward  and  have  been  studied  extensively.  Although  the  theory  may  predict  certain 
thermal  characteristics,  significant  uncertainties  still  exist  in  clinical  applications  resulting  in  less 
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than  desirable  correlation  between  onscreen  TI  values  and  actual  temperatures  within  the 
biological  medium.  However,  TIs  were  designed  to  be  conservative  and,  for  the  informed  user, 
onscreen  TI  values  can  be  useful  in  determining  the  potential  for  heat  damage  during 
examinations. 

Wavefront  “Movies” 

From  a  display  viewpoint,  X-Z  multiframe  “movies”  can  help  in  understanding  the  qualitative 
nature  of  ultrasound  wave  interaction  with  the  medium,  objects,  and  other  structures.  Such 
models  exist  in  literature  and  would  make  a  useful  addition  to  NUPROP.  For  instance,  Wojcik, 
Mould,  Ayter,  and  Carcione  [1998]  use  a  model  based  on  the  pseudospectral  code  that  eliminate 
the  use  of  lower  order  space  and  time  derivative  approximations  which  ultimately  lead  to  greater 
versatility  and  efiSciency  at  the  expense  of  numerical  accuracy.  Also,  their  model  relies  on 
accurate  modeling  of  abdominal  wall  morphology  to  capture  the  critical  phase  aberrant  nature  of 
biological  interfaces.  Such  data  exists.  Careful  measurements  of  abdominal  wall  section  was 
accomplished  by  Hinkelman,  Liu,  Metlay,  and  Waag  [1994]  and  Hinkelman,  Mast,  Orr,  and  Waag 
[1997]. 

The  pseudospectral  method  provides  an  extremely  accurate  approximation  of  spatial 
derivatives  in  homogeneous  media  using  FFTs  on  a  uniform  Finite  Difference  (FD)  grid.  For  a 
more  detailed  treatise  on  the  pseudospectral  method,  see  Fomberg  [1996].  The  price  of  high 
spatial  accuracy  is  the  necessity  for  space  periodic  construct.  In  other  words,  solutions  exhibit  a 
“wraparound”  effect  at  the  boundaries,  whereby  a  wave  exiting  one  side  of  the  computational 
space  appears  on  the  other  side.  One  way  to  circumvent  wraparound  is  invoking  Berenger’s 
[1994]  bovmdary  condition,  which  forces  boundary  values  to  be  “periodically  small”  at  the 
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boundaries.  Domain  periodicity  possesses  another  challenge  in  that  waves  are  introduced  as  initial 
conditions  rather  than  boundary  conditions. 

Accurate  time-domain  calculations  use  explicit  time  integrators  such  as  the  4*  order  Runga- 
Kutta  and  the  4"*  order  Adams-Bashforth.  The  pseudospectral  technique  uses  multiple 
frequencies  to  accurately  model  power  law  frequency  dependence.  A  least  squares  fitting 
procedure  chooses  model  parameters  for  an  optimal  fir  over  a  specified  frequency  range. 

The  pseudospectral  technique  facilitates  the  generation  of  2-D  simulations  such  as  X-Z 
multiframe  imaging.  Such  methods  aUow  the  user  to  watch  the  progression  of  a  wavefront 
through  a  series  of  piecewise  homogeneous  layers  as  illustrated  in  Figure  7-4.  Snapshots  of  a 
pulse  propagation  through  a  piecewise  homogeneous  abdominal  wall  is  presented:  Wojcik,  et.  al. 
[1998]. 
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Figure  7-4.  Snapshots  of  wave  propagating  through  piecewise  homogeneous  layer  models  (Left:  without  fine  structures,  Right: 
with  fine  structures).  The  fine  structure  added  to  the  model  on  the  right  is  based  on  actual  histological  slices  of  abdominal  wall 
structures.  The  fine  structures  introduce  considerable  diffuse  scattering  relative  to  the  ideal  model. 

As  is  expected,  the  fine  structure  produces  considerable  diffuse  scattering  as  illustrated  in 
Figure  7-4:  right  As  the  wave  propagates  through  the  fine  structure,  multiple  reflections  and, 
consequently,  constructive  and  destructive  interference  on  the  reflection  signal  from  these 
piecewise  sections.  This  brings  up  another  fi^litful  area  of  improvement:  development  of  more 
complex,  more  realistic  phantoms  of  tissue  models. 

More  Complex  Phantoms 

However  accmate  a  propagation  model  is,  its  “Achilles”  heel  is  stiU  the  fidelity  of  the  tissue 
model.  More  realistic  models  are  crucial  for  simulating  ultrasound  images  for  qualitative  and 
quantitative  diagnostics. 


104 


Piece-wise  Homogeneous  (modeling  inhomoseneity) 


As  previously  mentioned,  Wojcik,  et.  al.  [1998]  provides  several  piecewise  homogeneous 
models  for  abdominal  wall  and  multi-tissue  structures  based  on  histological  slices  and  measured 
acoustic  signals  (Hinkelman,  et.  al.  [1997]).  Figure  7-3  shows  the  piecewise  homogeneous  model 
used  by  Wojcik,  et.  al.  [1998].  Each  section  or  structure  possesses  its  own  set  of  tissue 


parameters  as  given 


Tissue/ 

Material 

P  3 

[kg/m^] 

V 

[m/s] 

B/A 

Attenuation 

[dB/cm/MHz] 

B 

Water 

_ Ic — SI. - J - 

1000 

1500 

5.0 

0.002 

2.0 

Fat 

928 

1427 

10.0 

0.75 

1.0 

Connective 

1100 

1537  1 

7.87 

1.125  ^ 

1.0 

Muscle 

1041 

1571 

7.5 

0.55 

1.0 

Liver 

1050 

1577 

6.75 

0.4 

1.0 

Currently,  NUPROP  will  model  each  section  individually,  propagating  the  wavefront  from  the 
front  interface  to  the  back  interface  through  the  material  for  one-way  propagation  studies  (onset 
distance,  single  frequency  studies),  and  back  through  the  same  layer  for  two-way  propagation 
(image  formation  studies).  NUPROP  handles  “external  boundaries”  by  assuming  perfectly 
matched  layering  (PML).  Inter-tissue  boundaries  are  modeled  in  NUPROP  using  point  scatter 
difference  amplitudes.  In  other  words,  different  materials  possess  different  maximum  amplitudes, 
but  have  a  zero  mean  value  so  as  to  eliminate  artificially  enhanced  inter-tissue  boundaries. 


Phase  Aberrating  Interfaces 


Biological  tissues  have  the  potential  of  introducing  phase  aberrations,  or  temporal  shifts,  in 


ultrasonic  signals,  thereby  corrupting  the  reconstruction  of  the  image  (B-Scan).  One  method  for 
correcting  phase  aberrations  is  backpropagation  of  a  wavefront  given  known  experimentally 
measured  data,  employing  angular  spectrum  as  the  means  of  backpropagation.  Liu  and  Waag 
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[1997]  use  such  a  technique,  with  the  stipulation  that  the  inhomogeneous  medium  is  modeled  such 
as  by  a  phase  screen  placed  some  distance  away  from  the  aperture  in  an  effort  to  be  able  to  steer 
the  transmit  beam  effectively  to  other  regions  within  the  interrogation  space.  This  property  of 
beam  steering  is  crucial  in  image  formation.  Liu  and  Waag  [1994]  couple  the  backpropagation 
wave  with  a  reference  wave  to  generate  time-shift  compensation.  Clearly,  this  technique  is  highly 
dependent  on  the  accuracy  of  phase  measurement  or  phase  simulation.  NUPROP  currently 
assumes  no  phase  distortion. 
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Appendices 

Appendix  I.  Mathematical  Development  of  Bni^ers’ 
Equation 


Borrowing  from  J.  David  Logan  [1987],  given  the  Burgers’  relationship  in  (2)  and  substituting 
the  assumed  solution  form,  u  =f(z-vt),  letting  s  =  z-vt,  5delds 

dt  oz 

-  v/'(^)  +  fis)fXs)  -  yf\s)  =  0 

(I-l) 


Upon  integration,  (I-l)  becomes 
-  v/  (5)  +  is)  -  yf'is)  =  B 


(1-2) 


where  5  is  a  constant  of  integration,  noting  that  fis)f  is)  2 


=  _L[/(5)  -  (v  -  ^|v^+2B  )][f(s)  -  (V  +  Vv'+2B  )] 
ds  2v 


(1-3) 


Assuming  >  -2B,  employing  a  separation  of  variables,  and  integrating,  yields 

j  =  2r  J - ==J - ,  .  - - # 

[/(5)“(v-Vv^  +  2B)][/(s)  -  (v + \v  +2B)] 

_  1  (v  +  Vv^+2g)-/(x)  ^j.4^ 

2^v^+2B  f(s)-iv-ylv^+2B) 


for  (v-Vv^+2B)<f(s)<(v  +  Vv^+2g) 
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Solving  for  f(s)  yields 


/ -  ; -  -l/v^+2B 

(v  +  Vv^  +2g)  +  (v-Vv^  +2B)e’' 


l  +  e'’ 

(1-5) 


Thus,  the  traveling  wave  solution  to  Burgers’  equation  is  given  as 


u(z,t)  =  f{z-vt) 


i -  f— -  ^Vv^+2fi 

v  +  vv^  +2^ +[v-Vv^ +2^]g  _ 


— Vv^+2B 

l+e  ^ 

(1-6) 


To  avoid  confiision,  the  velocity  is  usually  denoted  as  c,  thus  (1-6)  becomes 


u{z,t)  =  f{z-ct) 

r- -  ri - 

c  +  Vc^  +2fi+[c-Vc^+2g]e 


LJL^c^+2B 

l  +  e  ■' 

(1-7) 


It  should  be  noted  that,  without  the  viscosity  term,  -vV^u,  the  wave  will  tend  to  “shock” 
faster  and  form  a  “breaking”  wave.  The  viscosity  term  acts  as  a  dampening  of  the  nonlinear  term. 


du 

u - . 

dz 
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Appendix  II.  Mathematical  Development  of  Fubini’s 
Solution 


Beyer  and  Letcher  [1969]  provided  the  following  derivation: 

The  particle  velocity,  m  =  ^,  is  a  function  of  density  (p)  or,  likewise,  condensation  (s)  such  that 

1 


p  =  p^{\  +  d^ldz)  or  5  = 


i\  +  d^ldz) 


-1 


(II-l) 


It  follows  that  M  can  be  written  as  a  general  function  of  d^jdz  such  that  u  =  Fid^jdz)  thus 


Such  that 


dzdt 


and  -^  =  /  ,  , 
dzdt  dz 


B  . 

Comparing  (II-3)  with (II-2),  letting  ^ 


/'  = 


± 


..A 

(i+a#/az) 


ai-2) 


(11-3) 


(11-4) 


Upon  integration  of  (II-4)  the  following  result  is  obtained 
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(\  +  d^!dz) 


B 

2A  ^ 


(II-5) 


In  the  absence  of  sound,  a^/dz  equals  0  and  m  equals  0,  such  that 


K  = 


±2Ac„ 


B 


(11-6) 


leading  to 


M  =  ±- 


2Ac^ 


B 


\ 


1-- 


(l  +  d^/dz) 


lAc  ( 

=  +±l£^  l- (1+5)2'' 
B 


J 


ai-7) 


By  restricting  the  wave  traveling  in  the  -i-z  direction  and  assuming  that  d^jdz  is  negligibly  small, 
then  the  rate  of  wave  propagation  £is  a  function  of  particle  velocity  becomes 


v  =  c. 


lA  , 

B 

1  + - 

V  2^4  y 


(II-8) 


The  general  solution  of  (6)  is  given  as 


u  =  F[(o{t  ±  — )] 

V 


(11-9) 


such  that 
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Equating  (11-13)  and  (II-l  1),  the  following  is  obtained 

1  d^u  _  d^u 
a>^  dz^ 


or 

d^u  _  2 

If  the  wave  possesses  a  boundary  condition  of 
m(0,0  =  w<,  sin(of) 

then  the  general  solution  will  be 


(11-10,11, 12,  and  13) 


(11-14) 


(11-15) 


(11-16) 
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M(z,0  =  MoSin 


0}\ 


V  y 


=  Un  Sin 


-u^  sin 


ojz  .  B  u 

cot - 1  + - 

2A 


x-(— +1) 


cot-kz\ 


1+ - 


2A  c 


o  J 


where  it  =  — .  By  using  binomial  expansion  on 


two  terms,  yields 
m(z,/)«w^  sin 


1  + - 


'lA  c. 


(11-17) 


and  keeping  only  the  first 


o  ) 


f,  f,  B]u]\ 

cot -kz 

l-h+ — 

1  1  2A\cJ\ 

(11-18) 

At  this  point,  the  concept  of  the  discontinuity  distance  (/)  must  be  introduced.  By 
observation,  one  can  see  that  a  discontinuity  will  occm  when  the  denominator  term  of  (II-5), 
d^jdz,  becomes  negative  1.  The  denominator  term  must  approach  negative  infinity  at  some 
propagation  distance  due  to  the  fact  that  large  values  of  u  will  propagate  faster  than  small  values 
of  u,  thus  giving  rise  to  wave  distortions.  This  phenomenon  is  commonly  referred  to  as  a  “shock 
fi-ont.”  Fubini  computed  this  distance  fi-om  the  origin  to  the  point  where  the  shock  forms  as 

1  = 


(Bl2A  +  \)Gm, 


(11-19) 


Substituting  (11-19)  into  (11-18)  yields 
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which  can  be  expanded  into  a  Fourier  series  such  that 


—  =  y  sin[«(fiY  -  kz)] 
K  »=i 


ai-21) 


where 


1  w 

B„=—  f — sm[n(a>t  -  kzy]d((ot  -  kz) 

'TT  J  17 


^  0  “o 


(11-22) 


Substitution  of  (11-20)  into  (11-22)  yields 


(11-23) 


Such  that 


Uo  7^1  rtz  \  I 


sin[n(<aif  -  kz)) 


(11-24) 

For  distances  in  excess  of  /,  this  explicit  solution  becomes  invalid.  Upon  examination  of  (11-24), 
several  characteristics  can  be  observed: 

1.  The  fundamental  component  of  the  wave  vvdll  decrease  as  the  wave  propagates,  as  well  as 


harmonics  increase 


2.  If  the  initial  wave  amplitude  is  significant,  then  the  wave  will  propagate  a  shorter  distance 
before  experiencing  shock 

3.  As  fi'equency  increases,  the  shock  distance  will  decrease 

4.  The  ratio,  t4/uo,  can  be  viewed  as  appropriately  weighted  sum  of  fi-equency  dependent, 
sinusoidal  basis  fimctions. 
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Appendix  III:  Linear  Frequency  Modulated  Input  Signal 
(Chirp) 

Another  multiple  frequency  input  signal  is  the  linear  frequency  modulated  input  signal 
commonly  referred  to  as  the  chirp.  Arguably  one  of  the  most  useful  linear  filtering  techmques  is 
the  use  of  the  quadratic  phase  filter  or  chirp.  The  primary  advantage  of  using  a  chirp  signal  over  a 
single  frequency  sinusoid,  or  apodized  sinusoid,  is  the  ability  to  have  a  very  wide  bandwidth, 
infinite  in  fact,  while  maintaining  reasonable  signal  gain  (or  signal  to  noise  ratio).  From  Easton 
[1997],  the  quadratic  phase  filter  has  an  impulse  response  of 


1 


(III-l) 


where  □  is  the  chirp  rate,  and  a  frequency  domain  form  of 


1  ±ijri-f 


(III-2) 


Hence,  the  frequency  domain  representation  of  the  chirp  is  a  chirp.  This  means  that  both  the 
impulse  response  and  the  quadratic  phase  filter  have  constant  magnitude  and  infinite  support.  For 
large  values  of  a,  the  space-domain  chirp  will  oscillate  slowly  whereas  its  frequency-domain 
counterpart  will  oscillate  quickly. 

Because  the  chirp  is  an  all-pass  filter,  the  complex  correlation  of  its  impulse  response  is  a  delta 


function  as  demonstrated  using  the  filter  theorem  in  the  frequency  domain 

o  h'  (x)}  =  3{/i(x)  ♦  h*  (-X)}  =  •  3{/i*  (-x)} 

=  =  l(^)  ^  3-'  {1(^)}  =  S{x) 


(III-3) 


This  is  a  veiy  important  characteristic  because,  if  the  return  signal  is  comprised  of  a  complex 
summation  of  chirp  signals,  i.e.  the  introduction  of  a  chirp  ultrasound  into  a  biological  medium,  it 
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can  be  analyzed  using  an  appropriate  “chirp  mask”  to  determine  signal  content,  penetration  depth, 
etc.  A  chirp  mask  can  be  comprised  of  a  series  of  linearly  shifted  chirp  signals,  representing 
different  penetration  depths.  Figure  III-l  illustrates  a  simple,  linear  “chirp  mask,”  512  by  512, 
shifted  by  one  pixel  in  the  x  direction  for  every  y-direction  pixel.  For  this  illustration,  the 
maximum  mask  extent  is  20  cm  in  both  the  x  and  y  directions. 

Line  ar  Chirp  Mask 

512by  512 

Y/X;  1.0 


0.00  0.00  0.10  0.15  0.20 

Propagation  Distance  (m) 

Figure  III-l.  Linear  Chirp  Mask.  The  linear  ehirp  mask  is  comprised  of  a  “stack”  of  a  given  chirp  appropriately  shifted  in  the 
mapped  direction  as  a  function  of  shift  in  the  propagation  direction.  For  this  example,  every  shift  in  the  propagation  direction 
is  matched,  one  for  one,  by  a  shift  in  the  mapped  direction 

A  linear  chirp  mask  is  one  that  does  not  introduce  nonlinear  effects  as  a  function  of 
propagation  distance.  So,  as  the  chirp  pulse  propagates,  no  harmonics  are  introduced. 
Additionally,  no  amplitude  modulation  is  introduced  by  the  medium  for  this  example.  A  return 
signal  from  various  regions  within  the  propagation  path  will  be  a  complex  summation  of  the  chirp 
signnf  appropriately  shifted  in  space  to  represent  the  distance  traveled  by  the  chirp  pulse.  The 
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complex  return  chirp  signal  is  correlated  with  the  chirp  mask  and  the  correlation  peaks  are 
computed.  For  this  example,  a  complex  return  signal  is  generated  from  the  summation  of  chirps 
at  positions  (chirp  numbers)  40,  70, 100, 130, 160, 190,  220,  and  250.  These  positional  numbers 
correspond  to  0.016,  0.027,  0.039,  0.051,  0.063,  0.074,  0.086,  and  0.098m  propagation 
distances,  respectively.  Figure  III -2  shows  the  effect  of  additive  noise  on  the  correlation  of  the 
return  chirp  signal  with  the  linear  chirp  mask. 
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Figure  III-2.  Effect  of  noise  on  correlation  between  complex  signal  and  chirp  mask.  The  top  graph  depicts  a  complex  return 
signal  with  no  additive  noise  and  its  correlation  as  a  fiinction  of  chirp  number.  Each  peak  matches  exactly  the  location  in  the 
propagation  path  of  each  reflection  point.  The  center  graph  depicts  a  noisy  return  signal  and  its  correlation  with  the  chirp 
mask.  Even  under  these  stressful  conditions,  the  peaks  are  discemable.  The  bottom  graph  depicts  an  extremely  stressful 
sensing  condition:  S/N  of  0.5.  The  correlation  peaks  are  essentially  obliterated  by  the  noise. 


For  no  noise  conditions,  the  correlation  peaks  are  clearly  distinguishable  and  correspond 
exactly  Avith  the  chirp  numbers  40,  70,  100,  130,  160,  190,  220,  and  250,  respectively.  As  noise 
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increases,  the  correlation  peaks  become  less  and  less  prominent,  until,  at  a  noise  level  tvvice  the 
signal  level,  the  correlation  peaks  are  nearly  lost,  save  a  few.  Note  that  under  these  conditions, 
boosting  the  gain  of  the  signal  channel  will  necessarily  boost  the  noise,  thus  providing  no  signal  to 
noise  improvement.  The  chirp  provides  high  gain  in  noise  or  clutter  environments  due  to  its 
autocorrelation  characteristics  (i.e.  impulse  response)  while  maintaining  an  infinite  bandwidth,  in 
theory.  In  practice,  the  chirp  still  provides  high  gain  characteristics,  but  the  apodization,  or 
spatio-temporal  windowing,  will  limit  the  bandwidth.  Likewise,  the  transducer’s  spatio-temporal 
characteristics  will  limit  the  actual  detected  (and  processed)  signal. 

Now,  add  the  complexity  of  nonlinear  or  harmonic  generation.  As  the  chirp  propagates  in  the 
medium,  each  fi’equency  within  the  chirp  duration  will  potentially  undergo  nonlinear  wave 
formation.  So,  the  question  arises:  can  the  Fubini  or  Burgers’  approaches  be  modified  to  account 
fi5r  a  linear  fi'equency  modulated  input  as  its  initial  condition?  Consider  the  fisllowing  derivation: 
Given  a  fiinction  f(z,t)  such  that 

(III"4) 


Then  the  second  partial  derivative  of f(z,t)  with  respect  to  z  is  given  as 

=  /  (z,  t)[-  ilTdc^  +  (ilTrkari  -  ilTdc^z)^  ] 
dz 


(111-5) 


And  the  second  partial  derivative  of  f(z,t)  with  respect  to  t  is  given  as 
=  /  (z,  t)[-  /2;k£)^  +  (ilnkoiz  -  ilnw^tY  ] 


(III-6) 


Now,  suppose 
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(III-7) 


de  ^  dz^ 


with  k  =  —  and  z  =  cj,  where  c„  is  the  group  wave  veloeity.  Then 

^  =  f(z,t{-i2jr(c„kf  +(i2m:lk^t-i2m;lkhf  1 

dr  ^ 

=  }f  (z,oj^“  +  {i2nc^k^t  -  i27!Cok^tJ  j 

=»  -i2m;lk^  =  y{-  i2;ik^) 


r  =  cl 


(III-8) 


in  the  absence  of  nonlinearities,  yielding 


dr 


dz 


(III-9) 


Comparing  (111-9)  with  (3-6)  implies  that  the  Fubini  solution  for  the  chirp  should  be  similar  to  the 
single  sinusoid  input,  with  the  basis  functions  now  being  chirps  with  varying  initial  frequencies 

such  that 

«  "  J  (nzlh 

—  =  2y  CHIRP{m)GA USS{ns) 

u„  h  (rtzll) 


where  CfflRP(ns)  is  defined  as  ,  GAUSS(ns)  as  ,  5  and  crdefined  as  previously. 

It  is  this  hypothesis  on  which  NUPROP  constructs  the  Chirp-nonlinear  result. 
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